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ABSTRACT 


This  report  derives  equations  predicting  satellite 
ephemeris  error  as  a  function  of  measurement  errors  of 
space-surveillance  sensors.  These  equations  lend  them¬ 
selves  to  rapid  computation  with  modest  computer  re¬ 
sources.  They  are  applicable  over  prediction  times  such 
that  measurement  errors,  rather  than  uncertainties  of 
atmospheric  drag  and  of  Earth  shape,  dominate  in  producing 
ephemeris  error.  This  report  describes  the  specialization 
of  these  equations  underlying  the  ANSER  computer  program, 
SEEM  (Satellite  Ephemeris  Error  Model).  The  intent  is 
that  this  report  be  of  utility  to  users  of  SEEM  for  in¬ 
terpretive  purposes,  and  to  computer  programmers  who 
may  need  a  mathematical  point  of  departure  for  limited 
generalization  of  SEEM. 
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I.  INTRODUCTION 


A.  General  Overview 

Earth-satellite  ephemeris  estimation  (i.e.,  position 
prediction)  is  fundamental  to  many  space-related  operations. 
Measurements  by  friendly  space-surveillance  sensors  are 
computer  processed  to  yield  necessary  ephemerides.  Each 
ephemeris  thus  provided  has  some  characteristic  accuracy. 

The  prediction  of  ephemeris  error  is  also  important, 
both  operationally  and  in  the  planning  of  space  surveil¬ 
lance  systems  and  of  data  reduction  procedures. 

The  problem  addressed  here  is  the  mathematical  predic¬ 
tion  of  ephemeris  error,  as  it  results  from  measurement 
error  alone.  The  results  are  valid  under  conditions  where 
one  may  validly  ignore  uncertainties  of  atmospheric  drag  and 
of  Earth  shape.  A  major  requirement  was  that  resulting  equa¬ 
tions  be  suitable  for  computer  programming  to  obtain  rapid 
calculations  with  modest  computer  resources. 

This  report  presents  a  detailed,  general  parametric  solu¬ 
tion  to  the  above  problem.  This  report  gives,  in  particular, 
the  somewhat  specialized  form  of  that  solution,  which  is  the 
basis  of  the  new  ANSER  computer  program,  SEEM  (Satellite 
Ephemeris  Error  Model).* 

SEEM  demonstrates,  with  a  time-shared  HIS-635  computer, 
the  requisite  programming  suitability  of  the  mathematical 
results  herein.  Reference  1  describes  successful  use  of 


*  A  FORTRAN  program,  written  by  the  author  of  this  report, 
as  yet  unpublished. 
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SEEM  empirically  to  investigate  conditions  of  its  validity 
for  ephemeris  error  prediction  vis-a-vis  Earth-based  radar 
sensors.* 

This  report  is  directed,  first,  to  users  of  SEEM  who 
may  wish  to  understand  its  fundamental  interpretation.  It 
is  also  directed  to  computer  programmers  who  may  wish  to 
modify  certain  specializations  of  the  current  version  of 
SEEM  and  need  a  mathematical  point  of  departure.  The  in¬ 
tent  is  that  material  here  be  accessible  to  engineers  and 
scientists  who  do  not  necessarily  specialize  in  either 
statistics  or  astrodynamics . 

Thus,  this  report  is  semi-tutorial  in  style,  and  is 
derivationally  self-contained  to  the  extent  practical.  Some 
original  mathematical  developments  are  included  and  appro¬ 
priately  noted. 

B.  Technical  Overview 

The  purpose  of  this  section  is  to  provide  not  only  an 
overview  of  report  organization,  but  also  a  substantive  dis¬ 
cussion  of  ephemeris  error  estimation  sufficient  (1)  for 


*  A  summary  of  established  SEEM  applicability  is  as  follows. 
The  model  validly  accommodates  drag  forces  for  satellite 
altitudes  above  about  180  km,  over  time  intervals— encom¬ 
passing  both  sensor  measurements  and  the  prediction  times 
—of  up  to  at  least  9  hours.  The  model  validly  accom¬ 
modates  noncentral-forces  gravitational  fields  for  low- 
altitude  satellite  passes  by  as  many  as  three  Earth-based 
radars,  over  somewhat  longer  measurement-and-prediction 
time  intervals.  Assumptions  here  are  (1)  current  capa¬ 
bility  for  predicting  drag  forces;  (2)  current  under¬ 
standing  of  geoid  and  other  gravitational  perturbations; 
and  (3)  no  radical  radar  accuracy  improvements  beyond 
today's  state-of-the-art. 
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routine  interpretation  of  SEEM  inputs  and  outputs,  and  (2) 
for  appreciation  of  options  for  limited  generalization  of 
SEEM. 

Three  subsections  follow.  The  first  gives  the  overall 
technical  approach.  The  second  overviews  report  organiza¬ 
tion,  and  in  particular  that  of  Chapter  II.  The  third 
discusses  substantively  Chapter  III,  the  heart  of  the  report. 

1 .  Approach  to  Solution 

The  overall  approach  to  solution  oi  the  stated  problem 
is  to  ignore  drag  and  geoid  uncertainties  by  estimating 
ephemeris  errors  under  the  Keplerian  (central-force-field) 
approximation.  The  rationale  for  this  approach  is  that  the 
distance  of  a  "Keplerian"  satellite  from  a  Kepler ian-estimated 
position  should  approximate  the  distance  of  a  "real-world" 
satellite  from  a  position  estimated  via  perturbation  theory, 
provided  that  calculational  treatment  of  perturbations  is 
exact. 

The  cited  investigations  of  Reference  1  appear  to  con¬ 
firm  the  validity  of  the  Keplerian  approximation  for  the 
problem  at  hand. 

The  specific  analytical  approach  of  this  report  is 
parametric,  utilizing  standard  linear-algebra  procedures 
in  a  covariance-matrix  formulation. 

2 .  Organization  of  the  Report 

This  report  constitutes  three  chapters,  plus  four  ap¬ 
pendices.  The  non-statistician  reader  may  wish  to  read 
Appendix  A  before  proceeding  to  Chapter  II.  Appendix  A  is 
a  purely  tutorial  review  of  covariance  matrix  theory. 
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Cnapter  II  reviews  the  linear-algebra  theory  by  which 
one  may  estimate  ephemerides — in  a  Keplenan  universe — 
from  sensor  measurements.  This  chapter  serves  also  to  de¬ 
velop  the  notation  used  later  on.  The  tneory  of  Chapter  II 
accommodates  a  variety  of  sensor  types,  sensor  basing  both 
terrestrial  and  on  satellites,  and  a  variety  of  "unbiased 
statistical  estimators." 

Section  II. A  defines  useful  coordinate  systems,  and  some 
matrix  transformations  among  finite  and  infinitesimal  (error) 
vectors  in  those  systems.  Section  II. B,  drawing  upon  Appen¬ 
dix  B,  deals  with  transition  matrices  between  astrodynamical 
state  vectors  and  between  error  vectors  associated  witn  those 
states . 

Section  II. C  derives  the  sensor  "observation  equation." 
Section  II. D  treats  the  iterative  differential  correction 
process,  which  transforms  observations  into  an  estimated 
state  vector  corresponding  to  an  arbitrary  epoch  (i.e.,  in¬ 
stant  in  time).  A  subset  of  the  components  of  this  vector 
constitutes  an  ephemeris  estimate. 

Section  II. E  discusses  error  minimization  criteria  and 
associated  unbiased  statistical  estimators.  Finally,  Section 
II. F  outlines  calculational  shortcuts  for  (1)  estimating 
state  vectors  for  many  epochs,  and  (2)  revising  state-vector 
estimates  so  as  to  exploit  newly  available  measurement  data. 

The  next  subsection  describes  Chapter  III  Doth  organiza¬ 
tionally  and  substantively,  and  also  defines  the  supporting 
role  of  Appendix  C. 
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3.  Ephemeris  Error  Analysis 

Chapter  III  deals  with  prediction  of  errors  in  ephemer- 
ides  that  would  be  arrived  at  by  the  methods  of  Chapter  II.* 
This  error  analysis  assumes  perfect  convergence  of  the  iter¬ 
ative  differential  correction  process  of  ephemeris  estimation. 
Thus,  the  limiting  accuracy  characteristics  of  a  particular 
ephemeris  estimation  process  are  amenable  to  assessment,  even 
without  analyzing  the  practical  convergence  properties  of  tnat 
process. 

a.  Inputs,  Outputs,  and  General  Solution 

Section  III. A  defines  in  mathematical  terms  the  as¬ 
sumed  inputs  and  desired  outputs  of  the  problem. 

Key  desired  outputs  are  the  standard  deviations  of 
the  components  of  satellite  position  error,  for  the  time  of 
each  ephemeris.  These  are  to  be  expressed  in  a  coordinate 
system  selected  to  make  correlations  among  ephemeris  error 
components  vanish.  Further  desired  outputs  are  tne  orien¬ 
tation  angles  of  that  coordinate  system,  relative  to  a  "UVW" 
system  defined  as  having  the  following  axes  at  an  instant 
in  time: 


RADIAL 

" ALONG-TRACK ” 


Up:  geocenter  toward  satellite 

A  third  axis  orthogonal  to  the  other 
two  axes,  approximately  along  the  sat¬ 
ellite  velocity  vector  (exact  for  cir¬ 
cular  orbits) 


*  Actual  estimation  of  ephemerides  is  unnecessary  to  esti¬ 
mate  their  errors. 


CROSS-TRACK  -  Normal  to  the  orbital  plane,  di¬ 
rected  along  the  satellite  angular 
momentum  vector. 

The  above  standard  deviations  and  orientation  angles 
have  a  particularly  simple  physical  interpretation  if,  in 
addition  to  the  input  assumptions  listed  below,  the  probabil¬ 
ity  distributions  of  the  measurement  errors  are  of  multivariate 
normal  (Gaussian)  form.  Tnen  (see  Appendix  A)  one  may  inter¬ 
pret  the  standard  deviations  as  the  principal  naif-axial  di¬ 
mensions  of  a  "10"  ellipsoidal  confidence  volume,  oriented 
along  the  axes  of  the  rotated  coordinate  system.  One  may 
interpret  the  ellipsoid  as  centered  either  at  the  true  ephem- 
eris,  expressive  of  a  level  of  confidence  that  the  estimated 
ephemeris  will  fall  within  the  ellipsoid;  or  as  centered  on 
the  estimated  ephemeris,  expressive  of  a  level  of  confidence 
that  the  true  ephemeris  will  fall  within  the  ellipsoid. 

The  level  of  confidence  of  a  lo  ellipsoid  of  satel¬ 
lite  position  is  about  20  percent.  For  many  interpreta¬ 
tive  purposes,  a  3o  ellipsoid  (i.e.,  having  triple  the  dimen¬ 
sions  of  the  lo  ellipsoid)  is  more  useful,  providing  a  con¬ 
fidence  level  of  61  percent. 

Assumed  inputs  are  (1)  the  "true"  orbit  parameters 
of  the  satellite;  (2)  sensor  types  (measuring  any  subset 
of  the  quantities:  range,  two  angles,  and  their  respec¬ 
tive  rates)  and  locations  (sensors  may  be  fixed  or  may  move 
on  or  above  the  Earth's  surface);  (3)  sensor  envelopes  of 
geometric  coverage  (range  and  angle  extrema);  (4)  the  as¬ 
sumption  that  sensor  measurements  are  "unbiased";  (5)  "true" 
statistical  parameters  of  measurement  error,  in  the  form  of 


a  covariance  matrix  D  covering  all  sensor  measurements  made; 
(6)  the  covariance  matrix  C — some  approximation  to  D — 
characterizing  the  "estimator"  via  which  actual  estimation 
of  ephemerides  would  proceed;  and  (7)  the  instants  in  time 
corresponding  to  the  presumed  ephemerides  whose  accuracy  is 
in  question. 

There  are  still  further  assumed  inputs,  which,  how¬ 
ever,  derive  calculationally  from  above  inputs  (1),  (2),  and 
(3).  These  further  inputs  are  "ideal"  (error-free)  sensor 
observation  data  of  the  target  satellite.  A  "driver"  program 
to  generate  these  data  preexisted  SEEM  at  ANSER.*  This  re¬ 
port  does  not  give  the  full  mathematical  basis  of  such  a 
driver  program,  although  Sections  II. A  and  II. B  do  provide 
many  of  the  necessary  transformations. 

Regarding  inputs  (5)  and  (6)  above,  the  matrices  D 
and  C  typically  contain  >106  elements  each.  Apparently, 
their  general  parametric  specification  poses  a  practical 
difficulty.  Actually,  in  the  case  of  C,  this  difficulty 
must  in  some  sense  be  resolvable  relative  to  any  practical 
procedure  for  ephemeris  estimation,  since  any  such  proce¬ 
dure  involves  specification  of  C  (see  ensuing  discussion 
of  Section  III.C). 

Section  III.B  gives  the  general  solution  equations 
for  Keplerian  ephemeris  error  prediction.  These  equations 


*  A  multipurpose  Keplerian  program,  thus  far  unpublished. 
It  utilizes  a  modified  Earth  rotation  rate,  thereby 
correcting  to  first  order  for  the  drift  of  the  satellite 
orbital  plane  due  to  the  Earth's  equatorial  bulge. 
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provide  desired  outputs  as  a  function  of  assumed  inputs. 

These  equations  appear  to  present  a  further  practical  dif¬ 
ficulty,  in  that  they  require  extraction  of  the  inverse  of 
the  large  matrix  C." 

As  before,  this  difficulty  must  actually  be  resolv¬ 
able  for  error  prediction  relative  to  any  practical  proce¬ 
dure  for  ephemeris  estimation,  since  extraction  of  C_1  is 
necessary  there  also.  However,  this  difficulty  is  not  nec¬ 
essarily  resolvable  relative  to  "ideal"  ephemeris  estimation, 
for  which  one  must  take  C  *  D. 

b.  Detailing  of  the  General  Solution 

Section  II1.C  presents  candidate  representations  of 
D  that  resolve  the  specif icational ,  and  in  one  case  also  the 
inversion,  difficulties  identified  above.  These  represen¬ 
tations  serve  as  a  point  of  departure  for  selection  of  C 
representations  applicable  to  ephemeris  estimation  and,  con¬ 
comitantly,  to  ephemeris  error  prediction.  The  two  subsec¬ 
tions  of  Section  III.C  warrant  detailed  discussion  here. 

Subsection  III.C. 1  begins  by  assuming  that  in  the 
"real"  Keplerian  world,  raw  measurement  data  are  prepro¬ 
cessed  at  each  sensor  for  each  pass  as  follows,  for  entry 
into  the  ephemeris  estimation  process. 

Repeatedly,  raw  data  accumulated  over  an  interval  of 
a  few  seconds  are  suitably  averaged,  such  that  the  averaged 
results  correspond  to  the  instant  at  the  center  of  the  ob¬ 
servation  interval.  Known  corrections  for  systematic  error 


*  Practical,  not  theoretical,  invertibility  of  C  is  at  issue 
here.  Theoretical  existence  of  C“ ^  follows  from  the  defi¬ 
nition  of  C  as  a  covariance  matrix,  with  the  stipulation 
that  "perfect"  sensors  (i.e.,  having  any  sero  standard 
deviation  of  measurement  error)  are  not  allowed. 
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(e.g.,  sensor  calibration  corrections)  are  then  applied  to 
the  averaged  data  to  form  a  single  "observation  vector."  If 
the  sensor  happens  to  be  a  doppler  radar,  for  example,  ob¬ 
servation  components  would  be  range,  two  angles,  and  range 
rate.  At  the  end  of  the  pass,  the  collection  of  all  obser¬ 
vation  vectors  is  then  fed  into  the  ephemeris  estimation 
process.  Subsection  III.C.1  defines  D  to  be  the  "true" 
covariance  matrix  of  the  errors  of  all  observation  vec¬ 
tors,  over  all  passes  and  sensors. 

This  subsection  then  treats  the  errors  of  each  ob¬ 
servation  vector  as  the  sum  of  "noise  errors"  and  "residual 
bias  errors,"  the  latter  accounting  for  all  residual  system¬ 
atic  errors  in  the  observation.  By  assumption,  noise  errors 
may  be  correlated  with  each  other  within  an  observation,  but 
not  from  observation  to  observation.  By  further  assumption, 
the  noise  errors  are  uncorrelated  with  residual  bias  errors, 
within  an  observation  and  from  observation  to  observation. 

In  order  to  exclude  "perfect"  sensors,  all  noise-error  stand¬ 
ard  deviations  must  be  nonvanishing,  however  small.  In  order 
to  ensure  "unbiased"  measurements,  it  is  sufficient  to  assume 
that  both  noise-error  and  residual-bias  error  probability 
distributions  are  symmetric  about  zero. 

The  effect  of  this  decomposition  upon  D  is  to  render 
it  the  sum  of  a  "noise  matrix"  and  a  "residual  bias  matrix." 
Of  these,  the  noise  matrix  is  block  diagonal,  each  block 
being  the  covariance  matrix  of  a  single  observation  and  of 
dimensionality  at  most  6x6.  If  the  noise  errors  of  an 
observation  happen  to  be  uncorrelated  among  themselves,  the 
corresponding  block  will  be  diagonal. 
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The  residual  bias  matrix  may,  however,  be  relatively 
complicated,  with  widespread  off-diagonal  terms  representing 
long-term  correlations  among  residual-bias  errors.  Subsec¬ 
tion  II1.C.1  takes  a  first  step  toward  simplifying  this  ma¬ 
trix  by  assuming  zero  correlation  among  residual-bias  errors 
of  different  sensors.  Thus,  with  appropriate  organization 
of  D,  the  residual-bias  matrix  becomes  block  diagonal,  each 
block  corresponding  to  all  the  passes  by  a  particular 
sensor. 

Subsection  II1.C.1  concludes  by  developing  a  detailed 
parametric  representation  for  the  noise  and  residual-bias 
matrices,  structured  as  just  described.  Parameters  comprise 
various  error  standard  deviations  and  correlation  coefficien¬ 
cies,  with  general  functional  dependencies  upon  satellite 
position  relative  to  the  sensor. 

Thus,  Subsection  III.C.l  provides  D  structures  that 
are  physically  realistic  for  a  wide  variety  of  sensing  con¬ 
ditions.  It  also  provides  a  parametric  formalism  lending 
itself  to  practical  specification  of  D  as  an  input  to  ephem- 
eris  error  prediction.  However,  because  of  the  large  blocks 
of  elements  within  the  residual  bias  matrix,  the  structuring 
of  Subsection  III.C.l  is  not  generally  sufficient  to  provide 
practical  invertibilty  of  D.  Hence,  this  D-structure  does 
not  generally  permit  "ideal"  ephemeris  estimation  with  C  ■  D. 

The  objective  of  Subsection  III.C.2  is  to  furtner 
structure  the  residual-bias  matrix  so  as  to  arrive  at  an 
easily  invertible  D,  yet  accord  with  physical  reality  for  at 
least  some  measurement  circumstances. 
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Subsection  III.C.2  assumes,  for  any  given  sensor, 
that  the  residual  bias  errors  do  not  change  appreciably 
over  a  pass,  or  alternatively  over  several  passes  closely 
spaced  in  time  (a  "pass  mul t iplet" ) . *  This  subsection 
further  assumes  that  residual  biases  change  significantly 
between  pass  multiplets  (but  their  standard  deviations  do 
not  change)  such  that  residual-bias  correlations  vanish 
between  multiplets. 

Thus,  each  sensor  block  of  the  residual  bias  matrix 
decomposes  into  a  set  of  small  blocks,  each  corresponding  to 
a  pass  multiplet  for  that  sensor.  Each  multiplet  block  con¬ 
stitutes  a  set  of  partitions — corresponding  to  individual 
observation  vectors — that  are  identical  over  the  entire 
block. 

With  the  aid  of  a  derivation  detailed  in  Appendix  C, 
Subsection  III.C.2  infers  and  then  proves  the  validity  of  a 
closed-form  equation  yielding  D*1  as  a  function  (1)  of  the 
partitions  of  the  residual  bias  matrix  and  (2)  of  the  in¬ 
verses  of  the  partitions  of  the  noise  matrix.  (As  mentioned 
earlier,  these  partitions  are  of  maximum  6x6  dimensionality.) 
Hence,  except  for  matrix  multiplications  involving  the  parti¬ 
tions  of  the  residual-bias  matrix,  this  equation  reduces  the 
complexity  of  extracting  D-1  (as  structured)  to  that  of 
calculating  the  inverse  of  the  noise  matrix  alone. + 


*  Note  the  implication  that  residual-bias  error  is  insensi¬ 
tive  to  satellite  position  relative  to  the  sensor.  This 
assumption  may  be  inappropriate,  for  example,  if  it  should 
happen  that  atmospheric-refraction  uncertainties  become 
large  at  angles  near  the  horizon. 

+  This  equation  may  be  unique  to  this  report.  However,  a 
literature  search  was  not  feasible  within  the  scope  of 
this  analysis  effort. 
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To  summarize/  Section  III.C  provides  various  candi¬ 
date  D-structures,  including  structures  intermediate  to  those 
just  described,  for  use  in  detailed  expansion  of  the  general 
error-prediction  equations  in  Section  I2I.B. 

Section  III.D  continues  first  by  introducing  two 
broad  classes  of  C-structures  as  approximations  to  D,  and 
some  gradations  among  them.  Section  III.D  then  explicitly 
details  the  general  solution  equations  for  seven  combinations 
of  D-structures  and  C-structures. 

The  first  class  of  C-structures  constitutes  those 
congruent  to  the  noise  matrix  of  D,  i.e.,  those  wnich  are 
block  diagonal,  elsewhere  with  zeroes  for  every  element 
representing  correlations  from  one  observation  to  another. 
These  structures  are  readily  invertible  and  give  rise  to  what 
is  sometimes  referred  to  as  "weighted-least-squares"  (WLS) 
ephemeris  estimation.  WLS  estimation  ignores  all  correla¬ 
tions  from  one  observation  to  another. 

"Simple"  WLS  estimation,  a  special  case,  in  addition 
ignores  all  correlations  among  measured  quantities  within  an 
observation,  i.e.,  it  uses  a  C-structure  that  is  strictly 
diagonal.  This  is  equivalent  to  the  classical  estimation 
method  of  Gauss,  and  produces  an  optimum  fit  of  the  estimated 
orbit  to  actual  sensor  measurements  (see  Section  II. E). 

The  second  class  of  C-structures  consists  of  those 
allowing  nonvanishing  elements  that  represent  correlations 
among  observations.  These  structures  generally  incur  practi¬ 
cal  difficulties  of  inversion,  and  their  use  is  not  ordinarily 


2 


attempted  in  practice.  The  special  case  when  actually  C  *  D 
gives  rise  to  "minimum  variance"  estimation  and  produces — if 
calculationally  feasible — an  optimum  fit  of  the  estimated 
orbit  to  the  true  orbit. 

Because  of  its  ready  invertibility ,  the  "pass-multiplet" 
D-structure  of  Subsection  III.C.2  offers  the  opportunity  for 
minimum  variance  estimation  when  (1)  that  structure  is  valid, 
and  (2)  parameter?  of  measurement  error  are  known  with  suf¬ 
ficient  accuracy  tnat  •„  becomes  D. 

Section  5  provides  detailed  expansions  of  the 
general  error  prediction  equation  for  C-structures  that  are 
congruent  to  each  of  the  D-structures  of  Section  III.C,  and 
in  addition  for  the  WLS  C-structure,  which  is  congruent  to 
the  noise-matrix  of  them  all.  All  C-structures  are  distinct 
from  the  D-structures,  however,  in  that  their  parametric 
values  may  be  different — their  parametric  sets  may  even 
be  different  for  the  same  congruence  constraint. 

As  the  various  expansions  reveal,  the  amount  of 
feasible  detailing  of  solutions  is  quite  limited,  except 
for  those  C-structures  whose  inverses  can  be  extracted 
analytically.  Those  are  the  WLS  C-structure  and  the  "pass 
multiplet"  C-structure.  These  two  expansion  cases  comprise 
the  point  of  departure  for  the  specializations  of  SEEM. 

c.  The  Mathematical  Basis  of  SEEM 

Section  III.E,  the  final  section  of  the  final  chap¬ 
ter  of  the  report,  deals  with  SEEM.  Of  the  two  subsections, 
III.E. 2  presents  and  discusses  numeric  examples  of  SEEM  out¬ 
puts.  That  subsection  requires  no  further  discussion  here. 
Subsection  III.E. 1  gives  the  analytical  specializations  of 
SEEM. 
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As  input,  SEEM  accommodates  only  radars,  and  specifi¬ 
cally  only  those  that  operate  in  "altazimuth"  coordinates: 
azimuth,  elevation,  and  range. 

One  may,  however,  input  a  telescope-type  sensor  by  a 
strategem,  i.e.  ,  by  assigning  a  very  large  value  to  the  range 
measurement  error.  Thereby,  one  assigns  a  low  statistical 
weight  to  range  measurements. 

One  may  also  (to  some  approximation)  input  other 
range-and-two-angle  coordinates,  e.g.,  angles  relative  to 
the  boresight  of  a  phased-array  radar,  by  (1)  in  tne  driver 
program,  converting  to  altazimuth  coordinates  for  the  "ideal" 
observation  calculations;  (2)  in  the  driver  program,  finding 
a  geometrical  coverage  volume  in  altazimuth  coordinates  that 
approximates  the  true  coverage  volume;  and  (3)  assigning  angle 
errors  to  their  nearest  geometric  angle  analogs  in  azimuth  and 
elevation. 

SEEM  allows  correlations  only  among  errors  of  a  given 
measurement  component— e.g. ,  range-elevation  error  corre¬ 
lations  are  not  allowed.  Thus,  the  observation  blocks  of 
the  D  noise  matrix  are  diagonal,  as  is  each  small  partition 
of  the  residual  bias  matrix. 

Via  the  following  additional  assumptions,  SEEM  allows 
specification  of  the  error  performance  of  each  radar  in  terms 
of  six  parameters,  the  first  three  being  the  (constant) 
residual-bias  standard  deviations.  The  remaining  three  para¬ 
meters  are  the  noise  errors,  which  are  functionally  dependent 
upon  satellite  position  in  the  radar  field  of  view  as  follows: 

(1)  The  standard  deviation  of  azimuthal  noise  error 
is  proportional  to  1/cos  h,  where  h  is  elevation 
angle  (accounting  for  increasing  indeterminacy 
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of  azimuth  measurements  at  elevation  angles 
approaching  the  zenith).  The  constant  of  pro¬ 
portionality  is  hence  the  azimuthal  standard 
deviation  at  0°  elevation. 

(2)  The  standard  deviation  of  elevation  noise  error 
is  constant,  not  a  function  of  azimuth,  eleva¬ 
tion,  or  range. 

(3)  The  standard  deviation  of  range  noise  error  is 
constant,  not  a  function  of  azimuth,  elevation, 
or  range. 

In  the  light  of  (3)  above,  the  present  version  of 
SEEM  may  be  inappropriate  for  high-altitude  satellites,  where 
maximum  range  is  set  by  radar  range  performance  rather  than 
by  horizon-limited  line-of-sight.  (SEEM  validation  did  not 
include  satellite  altitudes  above  approximately  1,000  km.) 
SEEM  defines  all  pass  multiplets  as  containing  just  one  pass. 
Thus,  the  residual  bias  matrix  of  D — and  hence  also  D  itself— 
is  block  diagonal  in  one-pass  blocks. 

SEEM  provides  two  choices  of  C  for  characterizing 
the  ephemeris  estimation  process.  In  the  "minimum  variance" 
choice,  C  =  D.  In  the  "least  squares"  choice,  C  is  set  equal 
to  the  noise  matrix  of  D.  (SEEM  validation  was  conducted 
only  for  the  least-squares  choice.) 

SEEM  ephemeris  error  component  standard  deviations 
are  specified  in  UVW  coordinates  only,  and  do  not  include 
correlation  coefficients  that  may  not  vanish  in  those  coor¬ 
dinates.  A  further  output,  the  standard  deviation  of  the 
resultant  error  vector  is  invariant  with  respect  to 
coordinate-system  selection.  Hence,  the  UVW  resultant 
error  is  correct  even  without  coordinate  rotation. 
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Empirical  results  with  SEEM  indicate  that  in  fact 
the  error  ellipsoid  does  align  itself  with  the  UVW  axes 
soon — in  prediction  time — after  the  most  recent  pass 
(see  Subsection  III.E).  This  alignment  is  due  primarily  to 
the  effect  of  period  uncertainty,  which  makes  the  along- 
track  error  ordinarily  large  compared  to  radial  and  cross¬ 
track  errors. 
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II.  EPHEMERIS  ESTIMATION 


This  section  describes  linear-algebra  methodology  for  es¬ 
timating  satellite  ephemerides  from  sensor  observations,  as¬ 
suming  a  Keplerian  (central-force  field)  universe.  There  is 
no  restriction  as  to  selection  of  a  particular  statistical 
estimator,  except  that  it  be  "unbiased." 

Most  of  the  material  here  occurs — in  one  form  or  another — 
in  References  2  and  3.  A  first-order  correction  term  to 
the  error  transition  matrix  [i.e.,  in  Equation  (26)] 

does  not  appear  in  Reference  2,  but  is  well  known  in  esti¬ 
mation  theory.  The  explicit  representation  of  the  partial 
derivatives  of  may  well  be  new  as  developed  in 
Appendix  B,  but  are  available  elsewhere  in  somewhat  dif¬ 
ferent  form  (see  p.  B-6,  including  footnote]. 

To  the  exten:  practical,  the  notation  here  follows  the 
Herrick  standards  (Reference  3,  Aetrodynarical  Terninoloau , 
notation  and  Veage  (Appendix) t  pp.  477-511].  The  major  ex¬ 
ception  here  is  that  lightface  uppercase  Roman  letters  repre¬ 
sent  various  matrices  rather  than  specialized  astrodynamical 
quantities. 

A.  Reference  Frames,  Coordinate  Systems,  and  Transformations 

Let  the  time  tc  be  the  (arbitrary)  initial  epoch  of  the 
analysis.  Referring  to  Figure  1,  define  a  righthanded 
Cartesian  reference  frame  with  positive  z-axis  through  the 
Earth's  North  Pole,  and  positive  x-axis  intersecting  the 
Greenwich  meridian  as  it  happens  to  lie  at  tQ. 


*  This  x-axis  choice  promotes  algebraic  simplicity.  Con¬ 
version  of  ensuing  equations  to  a  system  with  x-axis 
positive  toward  the  vernal  equinox  is  straightforward. 


At  time 


t  ,  t  -  to 

let  the  satellite  position  be 


( 1 ) 


and  a  sensor  position  be 


(One  should  not  consider  the  sensor  position  as  necessarily 
on  the  Earth's  surface,  although  it  is  so  depicted  in  Figure  1 
for  ease  of  geometric  interpretation.) 

Again  referring  to  Figure  1  ,  define  a  topocentric 
(sensor-centered)  reference  frame  as  righthanded  Cartesian, 
with  positive  z'-axis  toward  the  zenith  and  positive  x'-axis 
toward  the  South  point  of  the  compass.  In  this  frame  let 
t.he  satellite  position  be 

x" 

P  =  y' 

. z 

With  these  definitions, 


r  -  rT  +  V 


(5) 
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FIGURE  1 

RELATIONSHIP  OF  INERTIAL  (xyz)  TO 
TOPOCENTRIC  (x’y'z')  REFERENCE  FRAMES 

(North) 
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FIGURE  2  e 

ALTAZIMUTH  COORDINATES 
IN  THE  TOPOCENTRIC  REFERENCE  FRAMES 


where  A  is  a  rotation  matrix.  Table  1  gives  a  representation 
for  Ft  an<3  A  in  terms  of  sensor  longitude  A,  latitude  ,  and 
geocenter  distance  v?.*  These  may  be  time-varying  quantities. 


Suppose  now  that  one  regards  the  components  of  p  as  func¬ 
tions  of  an  arbitrary  set  of  curvilinear  coordinates  qi,  q2»  q3» 
These  will  subsequently  become  the  angle  and  range  coordinates 
characterizing  operation  of  a  given  sensor.  (Their  particular 
significance  may,  however,  differ  from  sensor  to  sensor.)  Let 


Differentiate  Equation  (5)  with  respect  to  time,  obtaining 


r  *  rT  +  Kp  +  A  Jq 


(7) 


where  J  is  the  Jacobian  matrix  for  p  as  a  function  of  q.  Table 
1  gives  representations  for  A  and  J. 

0 

As  an  illustrative  example,  consider  the  case  of  a  sensor 
fixed  at  some  position  on  the  Earth's  surface  and  designed  to 
operate  in  altazimuth  coordinates.  Following  the  notation 
of  Figure  2,  one  may  write 


q 


h 

P 
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*  * 


(8) 


*  The  representation  of  A  is  a  standard  result.  One  may  de¬ 
rive  it  by  taking  products  of  elementary  rotation  matrices, 
which  provide  first,  a  rotation  of  the  primed  reference 
frame  about  its  y'-axis  through  the  angle  -  (*/2  -♦ ) ;  and 
second,  a  rotation  about  the  (now)  z-axis  through  the 
angle  -  (X  4  t). 


I 


J 
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TABLE  1 

TRANSFORMATION  QUANTITIES 
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and  arrive  at  the  specialized  transformation  matrices  of 
Table  2. 


Returning  to  the  general  case,  introduce  the  composite 
vectors 

It!  ■ 

which  will  subsequently  be  useful.  Equations  (5)  and  (7) 
provide  a  functional  relationship  between  these  vectors — 
unfortunately,  a  nonlinear  relationship  since  p(q)  is  nonlinear. 

However,  one  can  show  that  a  fully  linear  relationship 
does  exist  between  the  differentials  of  these  vectors,  of  the 
form 


These  vectors  will  represent  errors  at  a  particular  instant, 
so  that  in  performing  differentiations  t  is  to  be  held  constant. 

To  find  Q,  first  differentiate  Equation  (5),  regarding 
both  sensor  position  and  the  rotation  matrix  A  as  "known" 
(i.e.,  error-free)  and  hence  as  constants.  One  obtains 


6r  «  A  J  6q 

Differentiating  Equation  (7), 


(10) 


<f  -  AJ6q  ♦  AK6q  +  A  J  iq 


(ID 
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defining  a  new  matrix  K  via  the  relation 

A  6  J 4  =  A  K  6q  .  (12) 

From  Equation  (1)  one  can  derive  the  representation  for  K 
given  in  Table  1 . 

For  the  example  of  the  Earth-based  altazimuth-type  sen¬ 
sor,  one  can  further  derive  the  specialized  representation 
of  K  given  in  Table  2. 

One  can  now  find  the  general  form  of  Q  by  comparing 
Equations  (9),  (10),  and  (11).  Table  1  gives  the  result 
and  also  the  form  of  0“^.  One  can  prove  the  correct¬ 
ness  of  Q"1  by  taking  the  product  QQ"1. 

One  further  transformation  will  prove  useful: 


Here  L  is  the  rotation  matrix  taking  the  inertial-frame  (xyz) 
representation  of  r  into  a  L’VW  representation  (x"y"z")  de¬ 
fined  as  having  the  unit  vector  U(x"-axis)  directed  radially 
outward  from  the  Earth's  center  toward  the  satellite;  the 
unit  vector  W( z"-raxis )  directed  along  the  angular  momentum 
vector,  normal  to  the  orbital  plane;  and  hence  the  unit 
vector  V(y"-axis)  approximately  along-track  in  the  direction 
r— exactly  along-track  for  circular  orbits. 

To  find  L  in  terms  of  the  inertial-frame  (xyz)  represen¬ 
tation  of  r,  begin  by  defining  a  quantity 

•  -  r  x  r  .  (14) 
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(vector  cross-product),  proportional  to  the  angular-momentum 
vector.  Clearly  in  the  inertial-frame  representation. 


U  = 


i 


W  = 


_1 

s 


; 


v 


TsxTT 


■  (*  X  r)  ■ 
(s  x  r)y 
X  r>z- 


(15) 


(16) 


(17) 


But  the  above  nine  unit-vector  components  are  just  the 
direction  cosines  among  the  axes  of  the  two  reference  frames, 
and  hence  are  the  elements  of  the  rotation  matrix  relating 
the  frames.  That  is. 


where  notationally  U*  is  the  adjoint  of  U,  etc.  Hence,  L 
may  be  evaluated  from  Equations  (12)  -  (15).  Since  L  is  a 
rotation  matrix, 


L 


(19) 


B.  State  Vectors  and  Related  Transformations 

One  may  represent  the  six  parameters  of  a  Keplerian  orbit 
(including  the  instantaneous  position  of  the  satellite  in 
that  orbit)  by  the  state  vector 


# 


(20) 
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comprising  the  partitions  r,  f.  The  subscript  denotes  time 
functionality.  For  a  satellite  in  Keplerian  orbit,  the 
matrix  transformation 


x 


j 


X. 

1 


exists.  The  subscripts  denote  not  matrix  elements,  but 
rather  the  times  tj  and  tj.  The  propagation  m atnix  4>jj 
has  the  form 


where  I  is  the  3*3  identity  matrix  and  Appendix  B  gives 
expressions  for  the  scalar  functions  f,  g,  f,  and  g. 

One  can  gain  some  appreciation  for  4>jj  from  the  special 
case  of  a  circular  orbit,  for  which  (dropping  the  explicit  I 
which  one  is  still  to  understand  as  present) 


J  Circular 


Orbit 


cos (Ej-Ei) 


I  I  A  A  ■ 

|  isxnCE.-E.) 

■  I - - - 


_-n  sin  (Ej-Ei)  j  cosfEj-E^ 


Here  n  is  the  satellite  angular  rate  [radians/unit  time), 
and  Ei,  Ej  are  eccentric-anomaly  changes  since  t  *  0. 

Important  properties  of  4>ji  are: 

Functional  Dependence 


(21  ) 


(22) 


(23) 


(24a) 
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Composition 


Inverse 


Determinant 


i*jt>  *  1  • 

Note  that  when  Ej  *  Ej,  ♦ji  reduces  to  the  identity  matrix. 

One  can  find  the  transformation  between  small  errors  in 
state  vectors  by  differentiating  Equation  (19)  (i.e.,  while 

A  A 

holding  times  tj[  and  tj  constant): 


6  * 


=  (+.  . 
31 


This  defines  the  important  propagation  trror  matrix  ^jj. 
Appendix  B  develops  an  explicit  representation  for  'J'jj  in 
terms  of  Xj  and  (Ej  -  £j). 

Table  3  gives  a  special  case  of  this  representation, 
which  in  full  generality  is  algebraically  lengthy.  This 
case  corresponds  to  a  circular  equatorial  orbit,  where  the 
satellite  happens  to  have  the  position  component  x  *  0  [see 

A 

Equation  (2)]  at  time  t*  (i.e.,  the  time  corresponding  to 

!i). 


(24b) 


(24c) 


( 24d  ) 


(25) 

(26) 
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TABLE  3 

THE  MATRIX  (SPECIAL  CASE! 


Important  properties  of  j  are: 
Functional  Dependence 


Composition 

(*kj  4  *k;j)  <*;ji  4  *3i>  '  *ki  4  *ki  > 
Inverse 


'V  4  *•»> 


-1 


*  ^  +  ♦  .  . 
*3  *3 


Determinant 
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i,  E.  =  E. 


(27) 


(28) 


(29) 


(30a) 


*  lUnbounded  as  |Ej-E^| 

increases  without  limit], 

Wh en  E j  m  ^ji  *  0, 

The  inverse  follows  from  the  composition  property  by 
setting  k  «  i  and  using  Equation  (24a). 


(30b) 


One  may  prove  the  composition  property  by  differentiating 
Equation  (24b)  as  it  operates  upon  x^: 


{ ‘♦ki*!1 


Now  apply  the  definitional  Equation  (26)  first  to  the 
lefthand  side#  and  then  repeatedly  expand  the  righthand 
side: 


(31  ) 
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<*ki 4 


*  4  (£  *Ji’  “i1  4 


*kjI(*ji  4  4  “*kJ>M 


4,  .  f  (4  .  .  +  4'..)  6*.]  4  4',  .  dx. 
k3 1  31  31  a  3 


*kjI(tji  4  '*'ji)s,i1  4  '*'kjI(*5i  4 


(*kj  4 '*'k3,(*ji  4 • 


(31  ) 


Since  5*^  is  arbitrary,  the  compositional  Equation  (2B) 
must  hold. 

C.  The  Observation  Equation 

* 

Assume  that  at  time  tj  a  sensor  measures  a  subset  of 
the  components  of  the  satellite  coordinate  vector  q  and  its 
coordinate  rate  vector  q — for  example,  azimuth,  eleva¬ 
tion,  range,  and  range  rate.  Define  (  as  the  vector  of  true 
values  of  the  measured  components  and  y  as  the  corresponding 
vector  of  measurement  results.  Then 

yi  "  4  *  <32) 

where  ifj  is  the  measurement  error  vector  and  the  subscripts 
note  the  time  of  measurement  tj. 
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One  it, ay  now  define  a  matrix  Mj,  characteristic  of  the 
particular  sensor  making  the  measurement  at  tj,  by  the 
relation 


(33) 


Then  one  may  write  Equation  (32)  in  a  more  general  form,  con¬ 
venient  for  further  development: 


*5 


♦  1. 


Table  4  gives  examples  of  Mj  for  various  types  of  sen¬ 
sors,  all  of  which  operate  in  altazimuth  coordinates  (see 
Equation  (8)1.  That  is,  the  time  tj  is  characterized  by 
sensor  type,  not  only  as  to  the  measured  component-subset  of 
<lj,  <ij  (specified  by  Mj),  but  also  by  the  curvilinear 
coordinates  that  qj  represents.  Not  all  sensors  operate 
in  altazimuth  coordinates,  of  course.  Despite  interpreta- 
tional  differences,  the  mathematical  <o*m  of  Equation  (34), 
and  thutz^ctd  dimiruionality  qj  therein,  is  reasonably 
general  no  matter  what  the  value  of  %j. 


(34) 


To  proceed,  assume  next  that  before  the  measurement 
there  exists  a  preliminary  orbit  determination  in  the  form 

of  a  state  vector  i^H).  Here  the  asterisk  denotes  an  esti¬ 
mated  value  and  the  parenthesized  superscript  denotes  an 
initial  estimate.  (Methods  for  preliminary  orbit  determi¬ 
nation  are  discussed,  for  example,  in  Reference  2,  Chapters 
12  end  13.) 
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TABLE  4 

REPRESENTATIVE  (£,  M)  PAIRS 


(35) 


Now 


,  given  one  may  calculate 


♦jj1’  = 


and  then,  using  Equation  (21),  XjC1).  Then  one  may  use 
Equat  ions  (5)  and  (7)  to  solve  for  estimates  of  qj  and  qj, 
Let  these  estimates  be  denoted  q**1)  and  qjC1).  Finally, 
form  the  estimate 


t*(l)  = 


=  M. 


(1) 


(1) 


(36) 


One  may  then  subtract  Equation  (36)  from  Equation  (34) 
to  obtain 


L,i 


u) 


K5 


W(in 


•  •*  (1) 
q-q  u' 


J3 


(37) 


where 


tyi U1  5  ’3  *  «d*U) 


(38) 


One  can  further  utilize  the  known  sensor  location,  together 
with  the  estimated  quantities  of  the  preceding  paragraph, 
to  estimate  the  matrix  (Qjd))”1  (see  the  formulas  of 
Table  1).  In  the  light  of  Equation  (9),  one  may  write  the 
definition 


*  (o*(1>) 


-1 


q  - 
q 


-  ,*<17 


Note  that  according  * o  Equation  (9), 
A*jU>  *  "  *j  * 


(39) 


(40) 
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correct  to  first  order.  That  is,  correction  terms  on  the 

* 

righthand  side  of  higher  order  in  the  difference  Xj  -  xj 
may  exist.  Using  Equation  (39),  Equation  (37)  now  becomes 


M. 


(1) 


*  *j 


HD 


Suppose  one  wishes  to  find,  eventually,  an  estimate  of 

A 

the  state  vector  at  an  arbitrary  time  t*,  an  estimate  that 
is  to  be  an  improvement  over  an  estimate  that  is  simply 


,M1)  .  *M1)XM1> 

i  *ik  k 

i 

One  then  may  expect  it  to  be  advantageous  to  introduce  a%i 
as  defined  by 


(42) 


**3 


(1)  . 


= 

'  J1 


(1) 


+  * 


(1) 


Di 


A x{1} 


(43) 


This  is  a  first-order  version  of  Equation  (26).  The  paren- 

(D* 

thesized  quantities  may  be  computed  from  xx  ,  derived  from 
Equation  (42).  The  quantity  axx  is  of  course  unknown,  from 
the  viewpoint  of  the  satellite  observer.  Its  estimation  will 
be  appropriate  subsequently. 

Using  Equation  (43),  one  finally  may  obtain  from  Equation 
(41)  the  observation  equation 


Here 


(44) 


(45) 
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ft> 


In  the  observation  equation,  note  in  summary  that  one 
computes  ayj^and  Tjj/^from  the  observation  vector  yj 
and  from  the  initial  estimate  x^  already  assumed  to  be 
available.  The  remaining  quantities  are  unknown,  although 
subsequent  analysis  will  assume  knowledge  of  certain  statis¬ 
tical  properties  of  »jj. 

Further,  note  that  for  another  observation  at,  say,  time 
m,  one  will  have  redefined  the  quantity  such  that  it  may 

have  second-order  differences  from  the  axj  of  Equation  (44). 
The  next  subsection  will  ignore  such  differences  in  developing 
an  iterative  procedure  which — if  convergent — will 
eliminate  their  impact  upon  an  ultimate  estimate  of  xj. 

D.  Iterative  Differential  Correction 

Define  the  following  composite  quantities  for  a  set  of 
n  observations: 
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The  composite  form  of  Equation  (44),  j  *  1,  2,  n,  is 

then 


♦  n 


Note  that  this  equation  places  no  restrictions  upon  the 
time  separation  of  the  sensor  observations;  upon  the  sensor 

A  A  A  A 

"mix";  or  even  that  necessarily  tj  <  t2  <  •*•  <  tn-i  <  tn. 

Suppose  now  that  one  can  find  an  utimate*.  matrix  W*^)  (i.e., 
a  function  of  x * ) * 1  * ,  which  yields  an  estimate  of  j^M)  in 
Equation  (49): 


and  which  obeys  the  constraint  (not  an  approximation 


M*(l) -Ml) 

i  Ti 


(The  righthand  side  here  is  an  identity  matrix.)  The  next 
subsection  will  provide  a  class  of  such  estimators. 


Whatever  the  specific  version  of  one  can  interpret 

the  result  of  Equation  (50)  as 


K{1,]‘  ■  [v»;u,]‘  . 


'  using  Equation  (40).  That  is  U*i)*  iE  approximately  an  esti- 

mate  by  which  the  originally  given  orbit  estimate  *<(1)  (de- 
£  * 
rived  from  via  Equation  (42) J  was  in  error.  One  might 

hope  that  an  improved  orbit  estimate  would  be 


■;<”  - .;«« +  [*•!««]*  . 
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One  can  now  repeat  the  preceding  process 
by  and  obtaining 


replacing 


(2)is(2) 


(54) 


Continued  iteration  leads  to  a  sequence 


x 


Ml) 
i  ' 


,M2)  XM3) 
i  '  i  • 


which  may  converge,  i.e., 


Lim 

k-*~> 


0  , 


(55) 


depending  upon  the  form  of  W  and  upon  the  accuracy  of  the 
initial  estimate  x*^1). 

Suppose  after,  6ay,  k  iterations,  one  stops  the  iterative 
sequence,  when  there  remains  the  estimated  error 


Then  the  corresponding  form  of  Equation  (50)  is 


«<kl  -  w‘<k>4, 

* 

<*>4  uj 

*  Ax  P1*  + 
X 

w*(k)-n  . 

(57) 

using  Equations 

(49)  and  (51). 

The  ^iteration  analog  of 

Equation  (40)  is 

now 

a  ij- 

• 

*i  ' 

(58) 
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a  first-order  approximation  that  may  be  quite  good  if  sub¬ 
stantial  convergence  has  occurred.  Substituting  into 
Equation  (57)  and  rearranging, 


£(k) 


(59) 


This  is  a  result  of  fundamental  importance,  since 
it  specifies  the  error  in  x^*) — i.e.,  the  orbit  estimate 
ultimately  obtained  from  the  entire  process  and  that  con¬ 
tains  within  it,  as  a  partition,  the  ephemeris  estimate 

Assume  now,  and  henceforth,  that  the  convergence  has  been 

such  that  is  small  enough  to  be  ignored.  Then  for 

simplicity  one  may  denote  as  simply  xj  and  W*^) 

A  * 

as  just  W*  (i.e.,  showing  its  functional  dependence  upon  \j), 
to  obtain 


(*i"V  *  wi  *  »7  •  (60) 

Later  sections  will  analyze  Equation  (60). 

E.  Estimators 

The  purpose  here  is  to  define  and  discuss — but  not 
actually  to  derive — a  class  of  estimator  matrices  from 
which  one  may  select  a  particular  member  for  use  in  Equation 
(50).  The  approach  here  is  first  to  introduce  two  specific 
members  of  this  class,  and  then  to  generalize. 

Consider  a  relation  of  the  form 

(61 ) 


»  ■  Tfi  +  ( 


* 


in  which: 


o  ft  is  the  true  value  of  an  m-vector  whose  estimate  ft* 
one  wishes  to  obtain 

o  w  is  a  known  "measurement"  n-vector  (n  >  m) 

o  £  is  a  random  "error"  n-vector,  whose  value  is 
unknown  but  some  of  whose  statistical  properties 
are  known 

o  T  is  a  known  matrix,  not  a  function  of  ft. 

One  wants  to  obtain  the  estimate  ft*  via  an  estimator  matrix  W 
in  a  relation  of  the  form 

A*.  *  W»w  ,  (62) 

under  some  designated  optimization  criterion. 

Moreover,  one  desires  the  property  that  if  £  is  unbiased, 
then  the  estimation  error  (ft*  -  ft)  is  also  unbiased — i.e., 
one  desires  that  W  be  an  "unbiased  estimator." 

The  unbiased  estimator  property  translates  into  a  simple 
mathematical  constraint.  Substituting  Equation  (61)  into 
Equation  (62), 

fi*  «  V7fi  +  .  (63) 


Then  if  and  only  if 


WT  -  1  , 

(the  m  x  m  identity  matrix), 

V  -  t 


(64) 


(65) 
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and  if  {  is  unbiased  then  so  is  ( ft *  -  ft) .  Thus,  the 
txactne.**  ccn&tnaint  of  Equation  (64)  is  a  necessary  and 
sufficient  condition  that  W  be  an  unbiOLiid  . 

To  define  a  weighted  least  squares  criterion  for  esti¬ 
mation  of  p*  ,  begin  by  defining  the  quantity 


i.e.,  a  noise-free  measurement  vector  corresponding  to  ft*. 
Thus  if  ft*  is  nearly  equal  to  ft ,  then  **  becomes  nearly  an 
ideal  measurement.  One  might  reasonably  ask  that  V  be  chosen 
such  that  the  magnitude  of  w*  -»  be  minimal.  One  might  also 
ask  that  those  measurement  error  components  corresponding  to 
very  accurate  measurements  be  accorded  the  most  statistical 
weight.  That  is,  one  might  require  that  W  minimize 

*  2 


where  (o^)?  is  the  known  variance  of  the  ith  measurement. 

By  carrying  out  an  appropriate  minimization  procedure 
while  observing  the  exactness  constraint  (see  Reference  2, 
pp.  201-203),  one  can  obtain  the  result 

WLS  S  ^Weighted  Least  Squares 

.  (T^Tl-Vc^  . 

Here  by  definition  Cls  is  the  n  x  n  diagonal  matrix  with 
nonvanishing  elements 


<CTe> 


<Vi‘ 
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Several  features  of  this  result  are  significant.  First, 
it  clearly  satisfies  the  exactness  constraint.  Second,  the 
statistical  properties  of  £  that  must  be  known  are  its  com¬ 
ponent  variances  (i.e.,  knowledge  of  the  form  of  its  proba¬ 
bility  distribution  is  not  necessary).  Third,  the  form  of 
WLS  wakes  the  estimation  result  invariant  with  regard  to 
selection  of  m -component  dimensional  scale  (e.g.,  km  or  NM) . 

(Note  that  if  the  weights  are  arbitrarily  set  equal 

to  unity  as  in  "ordinary  least  squares"  estimation,  dimension- 
scale  invariance  no  longer  holds.)  Fourth,  the  inverse  of 
C^s  is  trivial  to  find,  so  that  numerical  evaluation  of 
Kls  is  straightforward  even  when  the  number  of  measurements 
is  large. 

One  may,  however,  adopt  a  different  estimation  criterion, 
and  arrive  at  a  somewhat  different  result.  Suppose  one 
decides  to  minimize,  not  the  measurement  residuals,  but  the 
state-vector  residuals — i.e.,  the  individual  component 
variances  of  the  estimation  error  (ft*  -  fj) .  For  minimum- 
variance  estimation,  one  is  to  minimize  each  of  the  quantities 

*  i*l,  2,  •••  n  , 

again  observing  the  exactness  constraint. 

If  one  carries  out  an  appropriate  minimization  procedure 
(see  Reference  2,  pp.  185-192)  assuming  now  that  the  "true* 

(  covariance  matrix  D  is  ?-nown,  one  can  obtain  the  result 

WMV  “  ^ Minimum  Variance  (69) 
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Significant  features  of  this  estimator  are  as  follows. 

First,  it  clearly  satisfies  the  exactness  constraint. 

Second,  the  statistical  properties  of  {  that  must  be 
known  are  its  entire  covariance  matrix  (not,  as  for 
merely  the  diagonal  elements  of  D).  Third,  W^y  provides  a 
result  that  is  properly  invariant  with  respect  to  dimen¬ 
sional  scale  changes.  Fourth,  the  inverse  of  D  may  not  be 
trivial  to  find  when  the  number  of  measurements  is  large,  so 
that  numerical  evaluation  of  W^y  may  not  be  straightforward. 

Fifth,  the  variance  o{  each  component  (p*  '  p)  indeed 
minimum  ( o r  Vfyy,  ol*  compared  to  the  \fj*  -  fi)  -component  var¬ 
iance*  oi  any  other  estimator  (including  Wls)*  However,  how 
can  one  obtain  D  with  assurance? 

In  fact,  D  will  not  be  exactly  known  in  practice,  but  may 
be  approximated  by  some  matrix  C  that  must  be  real,  sym metric, 
invertible,  and  have  positive  diagonal  elements.  Then  the 
practical  estimator  will  be 

«  -  (tV1!)'1  TV1  .  ,70) 

This  reduces  to  V?ls  if  C  •  C^g,  and  W^y  if  C  *  D,  but  in 
fact  represents  a  eta**  o(  e*timator*  where  C  is  selectable. 

Ease  of  calculation  and  estimation  accuracy  both  depend  upon 
selection  of  an  appropriate  C. 

One  may  show  (see  Reference  2,  p.  202)  that  the  esti¬ 
mator  of  Equation  (70)  results  from  minimization  of  the 
quadratic  form 

Tp)  ,  (71) 

subject  to  the  exactness  constraint.  This  generalized  opti¬ 
mization  criterion  is  known,  somewhat  confusingly,  as  the 
weighted  /cast  square*  criterion. 
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With  regard  to  the  preceding  subsection,  the  following 
correspondences  hold  for  the  kth  iteration: 

(72) 

(73) 

(74) 

(75) 

(76) 

Note  that  since  a  value  for  is  assumed  as  an  input 

to  each  iteration,  is  a  known  quantity  as 

assumed  in  the  minimization  of  the  quadratic  error  expression 
of  Equation  ( 71 ) . 

The  preceding  discussion  has  not  addressed  three  key 
questions.  Does  convergence  occur,  in  the  iterative  dif¬ 
ferential  correction  process,  for  an  arbitrary  selection 
of  C?  If,  for  a  given  C,  convergence  does  occur,  does  it 
necessarily  yield  a  unique  x*?  (That  is,  does  Equation 
(60)  have  more  than  one  solution?)  If  xt  is  unique  for 

* 

a  given  C-,  what  is  the  minimization  criterion  to  which  xj 
corresponds? 

In  fact,  convergence  may  or  may  not  occur  for  a  given 
selection  of  C.  Further  discussion  of  this  topic  is  beyond 
the  scope  of  this  paper. 


w 


n 

*  (k) 

*  (k) 
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One  may  show,  however,  (see  Reference  2,  pp.  437-440), 
that  if  convergence  does  occur,  the  limit  x^  is  unique  for 
that  particular  C.  Moreover,  the  resulting  minimizes 
the  quadratic  error  expression  of  Equation  (71),  wherein  T 
is  now  to  be  interpreted  as  T(xi). 

Note,  in  closing,  that  all  of  the  analysis  of  this  sub¬ 
section  presupposes  that  the  correct  functional  form  of  T 
is  known. 

F.  Calculations!  Strategies  for  Ephemeris  Estimation 

Suppose  one  must  obtain  ephemeris  predictions  for  a  se¬ 
quence  of  times  t*,  k  »  1,  2,  ... — that  is,  suppose  one 
must  obtain  a  number  of  estimates  x^  from  a  given  set  of 
measurements.  Suppose,  moreover,  that  during  the  time  in¬ 
terval  of  prediction,  additional  measurement  data  occasionally 
become  available.  One  then  desires  to  obtain  a  revised  set 
of  predictions  x £  in  near-real-time  with  the  arrival  of 
new  data. 

Calculational  efficiency  of  ephemerides  now  becomes  an 
issue:  is  it  necessary  repeatedly  to  carry  out  the  full 
iterative  differential  correction  process  for  each  x£? 

Calculational  strategies  do  exist  to  alleviate  this  prob¬ 
lem,  at  least  under  some  circumstances.  The  first  strategy 
simplifies  calculation  of  a  set  of  x£  from  a  given  set  of 
measurement  data.  The  procedure  is  first  to  find  one  state- 
vector,  say  x*,  via  iterative  differential  correction,  and 
tr.en  repeatedly  to  utilize  the  relation 

*Jt  *  *ki  *i' 

(see  Equation  (42}]. 


(77) 
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One  would  hope  to  obtain  in  this  manner — independent  of 
the  choice  of  i — the  same  set  x£  as  by  direct  use  of 
iterative  differential  correction  for  each  x£.  A  proof  of 
this  equivalence  will  conclude  this  subsection.  The  strategy 
of  Equation  (77)  is  of  quite  general  utility,  involving  no 
restrictions  as  to  the  nature  of  C  employed  in  W  [see 
Equation  ( 70 ) ) . 

Two  further  strategies,  the  "Bayesian  filter"  and  the 
"Kalman  filter  without  driving  noise,"  do  involve  such  re¬ 
strictions.  The  purpose  of  those  strategies  is  to  minimize 
the  recalculation  of  x\  for  Equation  (76)  when,  from  time 
to  time,  new  measurement  data  become  available.  The  key  to 
their  utility  is  treatment  of  each  batch  of  new  data  as 
having  errors  uncorrelated  with  errors  of  all  previous  data. 

Thus,  these  strategies  are  useful  for  particular,  block- 
diagonal  forms  of  C.  For  such  forms,  these  strategies  provide 
estimation  results  more  expeditiously  than,  but  identical  with, 
complete  "brute  force"  re-estimation  of  x*  from  old  and  new 
raw  data  for  a  given  C.  Reference  2,  Chapters  10  through  12, 
contains  a  detailed  discussion  of  the  Bayesian  filter  and  the 
Kalman  filter  without  driving  noise.  Their  further  discussion 
is  not  appropriate  here. 

* 

The  promised  equivalence  proof  will  demonstrate  that  x^ 
is  identical,  whether  arrived  at  by  iterative  differential 
correction  or  indirectly  via  Equation  (77).  That  is,  sup¬ 
pose  that  iterative  differential  correction  yields  directly 
*  *  ^  * 
at  tfc  a  value  x^,  and  directly  at  t*  the  value  xj,  whence 

a  value  xj<  obtains  via  Equation  (76).  The  problem  is  to 
prove  that 

**  .  <78> 

*k  Xk  9 
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Now  frorc  Equation  (60) , 


* 

*  -  I 


* 

wi 


and 


*  • 


» 


t  * 

where  the  W-arguments  are  respectively  xj  and  xk'*  Given 
Equation  (77),  then  by  Equation  (26)  to  first  order 


(79) 


(80) 


*k  “  *k  ~ki  ^xi  “  * 


(81  ) 


where 


_  *  *  ,  * 

ki  vki  *  ki 


(82 


Premultiplying  Equation  (79)  by  I  ki  an<3  combining  with 
Equation  ( 81 ) , 


*k  "  *k 


*  * 


-Jci  Wi  •  ” 


(S3) 


Now  from  Equation  (70)  and  the  correspondences  of  Equa- 

* 

tions  (72)  to  (76),  the  general  form  of  Wj  is 

Wi  *  \Ti  C  Ti)  Ti  c  A  .  (84) 

Using  this,  one  may  expand  the  lefthand  side  of  Equation  (84) 
as  follows: 
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*  * 

ki  wi - 


* ki  (Ti  c  Ti )  Ti  c 


=  (^.;<r  [Kn-KirR^ 

-  [(v^RK  •  S*J]  K-^JV1. 

The  last  step  involves  substitution  of  two  subsidiary  rela¬ 
tions.  The  first  is 

_*-l  _* 

“ki  *  “ik  ' 

a  restatement  of  Equation  (29).  The  6econd  is* 

=  (T‘-l\+ 

v  ki;  \  ki  / 

Now  from  Equations  (45)  and  (47), 


(85) 


(87) 


(87) 


*  *  * 

Ti  •  *ik  -  Tk 


Substituting  this  into  Equation  (85),  one  finally  has 


*  * 

:ki  Wi 


*  (T 


vvr1 


W, 


(88) 


(89) 


*  This  relation  holds  for  any  invertible  matrix  A.  Taking 
the  transpose  of  both  sides  of 


one  has 


A  A"1  -  1 


(A"1)tA+-  1 


/a-1 \ * 


Hence  (A"1)  is  the  inverse  of  A  ,  as  in 

(A“1)+  «  (A+)"1 
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using  the  definitional  Equation  (84).  Upon  substitution  into 
Equation  (83), 


* 


(90 


a  relation  identical  in  form  with  Equation  (BO).  But  as  stated 
earlier,  such  an  equation  has  a  unique  solution,  and  Equation 
(78)  therefore  must  hold.  This  completes  the  proof  as  required. 


III.  EPHEMERIS  ERROR  PREDICTION 

The  preceding  chapter  describes  linear-algebra  method¬ 
ology  for  estimating  Keplerian  satellite  ephemerides  from 
sensor  measurements.  Specific  estimation  techniques  within 
this  methodology  depend  upon  selection  of  a  matrix  C,  which 
is  some  approximation  to  D,  the  covariance  matrix  of  measure¬ 
ment  errors  [see  Equation  (70)  f f . ] . 

The  purpose  here  is  to  develop  equations  for  estimating 
errors  in  the  ephemerides  that  would  be  arrived  at  by  the 
methodology  of  Chapter  II. 

Of  the  five  ensuing  sections.  Section  III. A  defines  the 
error  estimation  problem  in  terms  of  inputs  and  outputs. 
Section  III.B  derives  equations  of  the  general  solution, 
expressed  in  terms  of  the  covariance  matrices  C  and  D. 

Section  III.D  introduces  specific  representations  for  C  and 
D.  Section  III.E  details  the  equations  of  the  general 
solution  in  terms  of  those  representations.  Finally,  Section 
III.F  defines  and  discusses  the  ephemeris  error  equations 
underlying  the  computer  model  SEEM. 

The  general  solution  of  Section  III.B  here  is  well  known 
[see  Equation  (16),  Reference  4],  A  contribution  of  this 
analysis  is  the  representational  discussion  of  Section  III.C, 
and  specifically  the  analytical  matrix-inverse  given  by 
Equations  (143)  and  (144).  When  used  in  the  estimator  W,  it 
affords  calculational  efficiency  plus  some  accuracy  improve¬ 
ment  over  the  conventional  least-squares  approach  to  emphem- 
eris  estimation. 

The  ensuing  analysis  assumes  reader  familiarity  with 
covariance  matrix  theory  as  reviewed  in  Appendix  A. 
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i.  The  Error  Estimation  Problem 

This  section  defines  the  ephemeris  error  analysis  problem, 
in  terms  of  assumed  inputs  and  required  outputs. 

1 .  Assumed  Inputs 

Assume  that  the  following  inputs  are  available: 

a.  The  epoch  set  t^ ,  k  *  1,2,...,  for  which  ephemeris 
errors  are  desired 

b.  A  state  vector  xc  corresponding  to  some  epoch  tD 
(affording  a  complete  specification  of  the  "true" 
orbit  of  the  satellite) 

c.  For  each  tj  at  which  a  measurement  is  made,  the  set 
of  quantities  [see  Sections  II. A  and  II. C] s 

Qj'  Mj'  Aj»  * j'  rTj 

(affording  a  description  of  "true"  observables, 
the  subset  of  these  actually  observed,  and  the  sen¬ 
sor  position) 

d.  For  each  sensor  the  functional  forms  of  Table  1 
necessary  to  evaluate  the  matrix  Oj1  from  Qj, 
(affording  subsequent  evaluation  of  Tjj(Xj)  [see 
Equation  (45)] 

e.  Relative  to  the  measurement  error  vector  jj  :  its 
"true"  covariance  matrix  D;  its  covariance  matrix 
representation  C  used  in  the  estimator  W;  the 
assumption  that  rj  is  unbiased,  i.e.,  that  rj  ■  0. 
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The  prior  generation  of  input  c.  from  appropriate  sensor 
characteristics  is  a  standard  space-surveillance  problem  not 
within  the  scope  of  this  report  (although  the  equations  of 
Sections  II. A  and  II. B  are  useful  in  solving  that  problem). 

The  assumption  of  e.  that  tj  is  unbiased  is  subject  to 
the  following  interpretation.  Suppose  each  observation 
error-vector  tjj  [a  partition  of  rjj:  see  Equation  (34)] 
is  what  remains  after  application  of  calibration,  atmospheric 
refraction,  and  other  known  bicu  corrections  to  the  raw  ob¬ 
servation  data.  The  error  Vj  then  comprises  "noise"  and 
"residual  bias"  contributions  (see  Section  XII. C).  If  the 

probability  distribution  of  each  of  these  is  symmetric  about 
_  *  — 
zero,  then  tj j  —  0  for  each  tj  and  hence  v  *  0. 

2.  Desired  Outputs 

Let  Sfc  denote  the  covariance  matrix  of  the  error  (xfc  -  *k) • 

Let  rSjc  be  the  upper  lefthand  3x3  partition  of  Sj<;  i.e., 
let  the  elements 

(rSk)  =  (sk)  »  n,  n  «  1,  2,  3  •  (91) 

Bin  sin 


Then  *Sk  is  the  covariance  matrix  of  ephemeris  error,  since 

X 5  *  K  *  ‘k’oJ-v} 


rk  -  rk 


•  * 

rk  -  rk 


* 

r.  -  r, 


r  k"  rk 


»  E 


~  V  (,k  -  '  i^j  (,k  -  \>  (f>  - f 
.  *  -  - 


-  'it*  (rk  "  r  k>  +  ('k  “  (i*  ~  *  k)+ 
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(92) 


Desired  simulation  outputs  are 


a.  The  matrices  rSjc,  corresponding  to  the  desired 

* 

epochs  tfc,  k  *  1,2,... 

b.  "Error  ellipsoid"  interpretation  parameters  (orien¬ 
tation  angles,  semi-major  axes)  for  each  rS)t,  where 
orientation  is  specified  relative  to  the  L'VW  reference 
frame  (see  section  I. A], 


Note  (see  Appendix  A)  that  the  interpretation  of  rSfc  in 
terms  of  an  error  ellipsoid  centered  at  is  legitimate 
only  when  the  probability  distribution  associated  with 
( rk  ~  rk)  is  normal  and  unbiased.  This  condition  is  met  whenever 
the  probability  distribution  associated  with  rj  is  normal  and 
unbiased,  since  by  Equation  (6'')  the  error  (x£  -  x^) 

[containing  ( r £  -  r^)  as  a  partition]  is  a  linear 
transformation  upon  rj . 


B.  General  Solution 


The  purpose  here  is  to  derive  general  equations  giving 
desired  simulation  outputs  as  a  function  of  assumed  simula¬ 
tion  inputs. 


Consider  Equation  (60),  written  for  an  arbitrary  epoch  tj: 


*i  "  *i  *  *  • 


(93) 


Subsection  II. F  has  established  the  formal  invariance  of  this 
first-order  approximation  (see  Equation  (58)]  under  the  epoch 
transformation  Equation  (26),  itself  induced  by  Equation  (21). 

The  theorem  represented  by  Equation  (A-16)  allows  one  im- 
Med lately  to  write 


S 


i 


(94) 
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The  lefthand  side  is  a  desired  simulation  output,  but  the 

* 

righthand  side  depends  on  the  quantity  *i — according  to 
Section  III. A,  not  an  available  input. 

But  via  the  Taylor  expansion  in  vector  form, 

*  * 

Wi  =  wtx^ 

‘  M(»i>  +  jj^rj  (*1  -  v  . 


One  can  see  that  the  order  of  approximation  of  Equation  (93) 
is  preserved  in  writing 


*i  *  *i  *  Wi  * 


The  Wj^  here  is  a  function  of  *i,  not 


VL  =  W(*i) 


(T  +  cT1  T^'1  T*  C’1 


(see  Section  II. E),  with  a  column  matrix  whose  general 
partition  is 

(Ti)3  -  MJ  °3l  23i 


(98a) 


(see  Equations  (45)  and  (82)],  where 


SH  *  *Ji  *  +Ji  • 
Evaluation  of  Equations  (98)  la  for 


(98b) 


*3  "  *3°  *o  * 
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Thus 


K’i  D  WV+  ; 


(100) 


and  since 


ski<*i  -  "i1  • 


(101  ) 


one  has 


C  «  -r  c  r  + 

5k  “ki  5i  ”ki 


(102) 


Equations  (100)  and  (102),  together  with  Equations  (91),  (97), 
and  (98),  give  the  desired  outputs  a.  as  a  function  of  the 
assumed  inputs  a.,  b. ,  c. ,  d. ,  and  e. 

It  remains  in  this  subsection  to  obtain  the  output  b. 
from  rSfc  expressed  thus  far  relative  to  the  inertial  (xyz)  frame. 

By  Equations  (13)  and  (A-1B), 


N. 


L  rS)l 


(103) 


Diagonal ization  of£rSfcjuvw*  if  carried  out  by  an  appropriate 
numerical  procedure,  yields  eigenvalues  (oi)^»  (02)^1  (0  3 )  ^ 
and  corresponding  normalized  eigenvectors,  which  one  may  denote 
as  * 2*  e3*  According  to  the  analysis  of  Appendix  A, 
the  semi-major  axes  of  the  "l-o"  ellipsoid  have  values 

0  2*  0  3* 
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Or.e  way  to  interpret  the  ellipsoid  orientation  is  as 
follows.  Pick  the  first  eigenvector  and  denote  the 
angles  between  and  U,  V,  W  as  ®12»  anc*  ®13»  respec¬ 

tively.  Then 

cos  0^  «=  U  , 
cos  012  =  e1+  v  / 
cos  e13  =  "  f 


1 

0 

-  - 

0 

u  « 

0 

0 

,  V  «= 

1 

0 

,  W  * 

0 

1 

These  are  the  direction  cosines  of  the  orientation  of  the 
a i-axis  of  the  error  ellipsoid.  Note  that  one  may 
arbitrarily  change  the  sign  of  ei  if  ease  of  interpreta¬ 
tion  of  the  angles  is  improved  thereby.  Such  a  change  does 
not  upset  the  normalization  of  *1  and  physically  means 
that  one  may  take  either  of  the  "oj-ends"  of  the 
ellipsoid  to  be  "positive." 

One  can  interpret  each  of  the  remaining  eigenvectors 
similarly,  completing  extraction  of  desired  outputs  b. 

C.  Measurement  Covariance  Matrix  Structures 

The  purpose  here  is  to  introduce  some  candidate  algebraic 
structures  for  the  "true"  covariance  matrix  D.  These,  perhaps 
with  still  further  approximations,  are  also  candidate  struc¬ 
tures  for  C. 


mgmmmm 


(104a) 

(104b) 

(104c) 

(105) 
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The  first  of  the  following  subsections  introduces  struc¬ 
tures  by  consideration  of  the  "noise"  and  "residual  bias" 
concepts  mentioned  earlier.  The  second  subsection  further 
details  this  structure  into  a  form  which,  although  of 
limited  generality,  does  permit  analytical  extraction  of 
the  matrix  inverse— of  importance  since  one  would  like 
to  use  W  with  C  =  D. 


1.  Basic  Structure 

Let  each  observation  epoch  t ^  now  be  denoted  as  “t j , 
where  a  indexes  the  sensor  making  the  measurement;  m  indexes 
satellite  "pass"  through  the  field  of  view  of  that  sensor; 
and  j  is  now  to  be  regarded  as  indexing  an  observation 
within  a  pass. 

Consider  the  structure  of  D  first  with  regard  to  parti- 

/•* 

tions  corresponding  to  each  observation  epoch  “t™,  and  then 
with  regard  to  the  "microstructure"  within  such  partitions. 


a.  The  General  Observation  Partition 

Let  the  observation  error  vector  arj^  correspond  to 
the  epoch  “t^.  Let  the  ordering  of  the  aD;j  within  the 
composite  vector  17  (see  Equation  (48))  be  hierarchic,  such 
that  a  varies  the  slowest,  m  the  next  slowest,  and  j  the 
fastest.  The  ordering  of  other  composite  vectors  and  of  the 
matrix  D  will  of  course  correspond.  Note  that  freedom  to 
select  this  indexing  hierarchy  exists,  since  up  to  this 
point  the  analysis  has  not  restricted  interpretation  of 

A 

the  sequence  tj  (see  II. D). 


Now  introduce  the  decomposition 

®l *  a*j  +  B<* 


t 


(106) 
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w’"  ■  the  requirement  that  the  term  °£  ™  account  for  any  cor¬ 
relations  that  may  exist  from  one  observation  to  another. 

n.  jri  Ct  rn 

That  is,  i>j  is  to  be  interpreted  as  a  "noise"  term  and  £j 
is  to  be  the  "residual  bias"  term  within  a^j'.  It  follows  that 

•*  *- v  • 

where 

Mathematically,  this  decomposition  entails  to  this 
point  no  loss  of  generality.  Physically,  the  desirable 
interpretation  is  that  the  %j1  result  from  random  receiver 
and  possible  random  external  noise  sources,  with  each 
"observation"  actually  deriving  from  some  small  data  set 
such  that,  for  appropriate  observation  spacing.  Equation 
(107)  holds.  This  physical  interpretation  clearly  implies 
some  practical  constraints  as  to  signal  processing,  and 
moreover  implies  that  always  aGjl  is  positive  definite. 

a  m 

Regarded  as  residual  bias  errors,  the  £j  by  contrast 
will  be  correlated  with  each  other  from  one  observation  to 
another,  at  least  for  observations  not  too  widely  separated 
in  time  and  made  by  the  same  sensor.  Thus,  one  might  set 


9 


assuming  no  correlations  from  one  sensor  to  another. 


(107: 

(108 


(109 

(  110 
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Then  by  the  symmetry  of  Equation  (110), 


ojiin  a„nm 
"jk  “kj 

Also,  will  be  non-negative  definite,  taking  into  ac¬ 

count  that  residual  bias  errors  may  sometimes  vanish. 

Finally,  assume  that  the  noise  errors  are  uncor¬ 
related  with  the  residual  bias  errors: 


(ill) 


‘IC'SMI  -  °  • 

(112) 


Taking  all  of  these  relations  into  account,  the 
general  observational  partition  of  D  is  then 

-  4„e  (“G”  6jk  *  “«™)  < 

with 

aB-mn  „  Ba_nm  \ 

V  *  Dkj  • 

The  resulting  D  is — as  required — real,  symmetric,  and 
invertible,  and  has  positive  terms  on  the  diagonal. 

Note  that  sufficient  conditions  for  tf  to  be  unbiased 

are 


I 


(113) 

(114) 


j 

(115) 


o  m 

•  'j 


-  0 


(116) 
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1 


and 


a  in 

£j 


(117 


Physically,  one  can  regard  these  as  averages  for  each  sensor 
over  large  ensembles  of  observations,  with  probability  dis¬ 
tributions  symmetric  about  zero. 

b.  Microstructure 

Now  consider  the  problem  of  representing  structure 
a  m  ,o  mn 

within  the  matrices  Gj  and  Hj^.  As  will  soon  appear, 
there  is  a  problem  in  establishing  a  reasonable  notation  in 
which  that  structure  may  be  specified.  This  problem  will 
receive  primary  attention  here.  Possible  functional  depend¬ 
encies  of  certain  quantities  will  receive  limited  considera¬ 
tion. 


oi  •  m  . 

and  q j ,  representing  a  range  and  two  angle  variables,  and  the 

rates  thereof.  Some  subset  (possibly  all)  of  these  six  quan¬ 
tities  is  actually  measured,  resulting  in  the  errors  of  Equa¬ 
tion  (106).  If  one  indexes  the  components  of  that  Equation 
by  p,  p  *  1,2,...,  £  6,  then* 


In  the  present  notation,  observables  at  “t™  are  aq™ 


(■•;),  ■  (*•;),  •  (••;), 


(ns: 


One  may  now  introduce  a  useful  representation  for 


u_iu  _  , , 

Gj  as  follows. 


*  The  index  p  relates  to  the  components  of  °q  j*  and  aq  j1  not 
directly,  but  via  the  sensor  characterization  matrix 
previously  denoted  Mj  (see  Equation  (33)]. 
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F 


a  in 


Define  a  matrix  of  standard  noise  deviations  v2j 
as  having  the  general  element 


(?"L  * 


where 


(>'),  •  [■IM.Ml" 


(119) 


(12D 


(Here  q  *  1,2,...,  an  index  having  nothing  to  do  with  the 
vector  q.) 


Define  a  matrix  of  noise  correlation  coefficients 

a  m 

vRj  as  having  the  general  element 


otherwise 


1/  p  =  q  ; 


pq 


(121a 


(121b 


Then  from  Equation  (108), 

“g?  .  °Z  w  •  aRin  .  a  T  w 
j  vj  Vj  V  *  j 

a  m 

This  representation  of  Gj  separates  the  standard  deviations 
of  the  components  of  °»< j'  from  the  correlations  among  them. 
One  may  expect  that  the  correlation  matrix  will  depend  upon 
the  control-system  design  of  the  specific  sensor  a.  xf  the 


(122) 
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am  am 

components  of  vj  are  uncorrelated,  reduces  to  the  iden¬ 

tity  matrix.  One  may  further  expect  that  the  size  of  the 
standard  deviations  of  error  components  will  depend  in  general 

upon  signal  strength  and  receiver  noise.  That  is,  one  expects 
u  in  ci  in  ex .  fti 

that  v^j  will  be  a  sensor-dependent  function  of  Qj,  **j, 

a 

and  perhaps  also  of  some  parameter  set  ft  associated  with  the 
target — e.g. ,  radar  cross-section. 

A  representation  for  “Hj'fc,  analogous  to  Equation 
(122),  is  obtained  as  follows. 

Define  a  matrix  of  residual-bias  standard  deviations 
m 


e  I  j  as  having  the  general  element 


Hq  - 


where 


(:•;),•[*  If' 0 


1/2 


Define  a  matrix  of  residual-bias  correlation  coef- 
a  mn 

ficients  cRj}<  as  having  the  general  element 

*  1  if  p  ■  q  and  m  -  n  and  j  *  k; 

VE  3h/pq 


0  if  (°°i)P  °r  (“°")p  * 


otherwise 


M 

P 

L-i 

EC 

(123) 


(124) 


(125) 


61 


Then  from  Equation  (11 1), 


°«S?  -  X  ■  >7  •  X  • 

As  before,  this  representation  separates  the  standard  devia- 

,  a  in  an, 

tions  of  the  components  of  £j  and  £  ^  from  correlations 
among  them.  If  the  standard  deviations  of  residual  bias  do 
not  change  with  time, 

am  an 

c“j  *  c2k  *  (127) 

However,  this  may  not  be  the  case,  as  when  atmospheric  re¬ 
fraction  corrections  are  imperfect  at  low  elevation  angles, 
am  am 

Then  e  1 j  is  a  function  of  qj,  and  Equation  (127)  is 
only  an  approximation. 

a  mn 

The  form  of  the  correlation  matrix  ERjk  will  depend 

upon  the  sensor.  Its  elements  will  tend  to  diminish  for 

|a~  m  a*.  n  i 

tj  -  tfc|»  For  small  time  sep¬ 
arations  one  may  expect  that 

o  Rmn  a  _mm 

cjk  K  c*jj  *  (128) 

a  m 

If  the  components  of  f-j  are  uncorrelated  among  themselves, 
a  mn  J  a  mm 

then  cRjk  will  be  diagonal  and  eRjj  will  be  the  identity 

matrix. 

2.  Further  Structure 

The  purpose  here  is  to  introduce  further  assumptions  about 
the  structure  of  D  that  promote  ease  of  inversion.  All,  some, 
or  none  of  these  further  assumptions  may  be  valid  for  a  given 
problem. 


€2 


Suppose  for  sensor  a  there  exists  some  time  interval 
a 

Tzero»  which  is  the  minimum  separation  between  observations 
over  which  residual-bias  error  correlations  vanish: 


a-irm 

tV 


(12s) 


Tnen  by  Equations  (114)  and  (126),  a  condition  on  the  general 
observation  partition  of  D  is 


o6D™  -  0,  |aij  -  nji  » 0 

(130) 

Of  interest  here  is  the  case  where  the  separation  be¬ 
tween  passes  by  sensor  a  occasionally  exceeds  “t  zero.  (That 
is,  the  assumption  is  not  made  that  every  n-pass  separation 
is  larger  than  aT2ero.) 

To  arrive  at  an  appropriate  mathematical  formulation 
of  such  a  situation,  define  an  index  L  =  1,2,...,  which 
counts  pass  separations  for  which 


a^m+1  _  o^m  ^  aT 


zero 


Here  atn'+1  is  the  first  observational  epoch  of  pass  m+1  for 
sensor  a,  and  at111  is  defined  to  be  the  last  observational 
epoch  of  pass  m. 


Such  pass  separations  decompose  the  a-pass  sequence 
into  a  sequence  of  "multiplets,"  each  comprising  one  or  more 
passes.  Hence  one  may  consider  L  to  be  the  index  of  de¬ 
coupled  pass  multiplets. 


(131) 
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Now  L  is  an  index  dependent  upon  o  and  m: 


L  ■  L (o»  o)  • 

One  may  express  the  functional  dependence  upon  m  recursively. 
Let 


L(a,  1)  =  1 


For  arbitrary  m,  if  Equation  (131)  holds 

L(a,  m  +  1)  «  L(a,  m)  +  1 


otherwise 


; 


L  (a  ,  m  ♦  1)  *  L  (a ,  in) 

Now  let 

aL„inn  -  /  a  .m  a,nt) 

"jlc  *  E  \  ei  'k  J  ' 

under  the  constraint  that  passes  m  and  n  are  in  the  same 
multiplet: 

L(a,  m)  »  L(a,  n) 


The  assumption  of  occasionally  decoupled  passes  is  just 

oHnn  .  ciL-Jim 

“jk  °L(a,  a),  L (a ,  n)  “jk 

(see  Equations  (126)  and  (129)]. 

Utilising  the  redundancy 


*  *L(a, 


b),  Kb,  n)  an 


(132a 


(132b 


(133a 


(133b 


(134a 


(134b 


(135) 


(136) 
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one  may  now  write  Equation  (114),  subsuming  Equation  (130),  as 

_  i  ,  /arm  ,  *  j  otLyTUiX 

°jk  6L(a,  m) ,  L(a,  n)  \  j  fimn  jk  jkj 


Thus  D  is  block  diagonal.  If  the  blocks  of  D  are  denoted  aLD, 
then  the  problem  of  inverting  D  reduces  to  that  of  inverting 
the  set  of  smaller  matrices  aLD. 


b.  Strongly  Coupled  Passes 


Suppose  for  each  sensor  a  there  exists  some  maximum 
time  interval aT  one,  for  which  residual-bias  error  standard 
deviations  and  correlations  are  invariant — i.e.,  Equations 
(127)  and  (128)  hold  when 


|°t“  -  at£| 


one 


In  particular,  suppose  °r  ©ne  is  greater  than  the 
duration  of  any  pass  multiplet  as  defined  in  the  previous 
subsection.  Then  Equation  (138)  holds  whenever 


(137) 


(138) 


L(a,  m)  «  L (a,  n) 


(139) 


and  over  each  multiplet  one  has  an  invariant  matrix 


oLu  =  aL„mn 
H  -  Hik 

One  may  refer  to  this  as  the  assumption  of  4tKongly  coupled 
paaaea. 


(140) 


A  necessary,  but  not  necessarily  sufficient,  condition 
for  the  existence  of  multiplets  that  are  each  strongly  coupled 
internally,  yet  are  decoupled  from  one  multiplet  to  another, 
is  that 


one 


« 


tero 


(141) 
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For  example,  if  at  zero  happens  to  be  one  satellite  period, 
then  each  multiplet  can  comprise  but  one  pass.  For  a  low- 
altitude  satellite,  if  “t  zero  ^  12  hours,  then  each  multi¬ 
plet  might  consist  of  a  pair  of  passes  on  successive  revolu¬ 
tions.  For  some  problems,  Equation  (141)  cannot  be  satisfied 
for  multiplets. 

When,  however,  Equation  (140)  is  valid,  substitution 
into  Equation  (137)  yields 


Djk  =  fia6  fiL(a,  m) ,  L(a,  6mn  6  jk  + 

When  D  has  this  structure,  the  general  observational  partition 
of  its  inverse  turns  out  to  be  available  analytically: 


(142) 


a& 


6oB  fiL(a,  m)  ,  L (a,  n) 


where 


These  expressions  for  a^(D_1)^jk  constitute  a  general¬ 
ization  of  the  results  of  Appendix  C.  That  appendix  derives 
the  inverse  of  a  matrix  having  its  general  partition  of  the 
same  form  as 


*Gj  6mn  6jk  +  ^ 


(143) 


(144) 
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[see  Equation  (142)],  but  with  the  restriction  that  °Gjl  and 
aIH  must  be  diagonal.  That  restriction  allows  those  two 
matrices  to  commute. 

However,  the  result  of  Equations  (143)  and  (144)  em¬ 
bodies  no  such  restriction,  as  one  may  show  by  proving  that 
those  Equations  yield  the  proper  result  for  an  inverse: 


r°s  •  ■*(■") 


‘aeWjk  • 


(145) 


The  proof  proceeds  by  direct  substitution  of  Equa¬ 
tions  (136),  (142),  and  (143)  into  the  lefthand  side  here: 

-X  V*?  4»*  V  •  *n  mb. 

+  Yti  S“Y  <L(°>  »>•  My.  *)  aL“  '  ^yB^k)  4in4ik 
"  Ji  Soy  4Ma,  My,  t)  °Lh  '  {yB  4My.  1),  MB,  nl^i)  <,LU  *L«f 


^mn  ^  jk 


-  4aB4Ma,»>,MB,n)  * 

*  4oB  SL(a,  K) ,  L(a,  n)‘ 01,11  '  (°®k) 

*  *•*(&  *•"<«*  l(“'  ’  SM«,*>.  Mo,  n)'(B°iJ  •  “Lu  •  0I,“)  • 

■HI1  • 


(146) 
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The  last  term  here  is  a  sum  over  just  the  multiplet  L,  and 

aL 

has  no  nonvanishing  terms  over  that  multiplet.  Since  B  and 
altJ  are  constants  over  the  multiplet,  one  can  factor  them  from 
under  the  summation.  Moreover,  one  can  see  (via  a  ’truth 
table"  analysis)  that 

L  (a,  m) ,  L(a,l)  L(o,  t) ,  L  (o ,  n)  L  (o  >  m)  L  (q  ,  n)  L  (o ,  m)  L  (q  ,  £■ ) 


Hence  one  may  continue  to  develop  Equation  (146)  as 
^mn  *  jk 

+  {o6  {L(a,ro)L(a,n)[-  °L“  +  + 


^a8  ^mn  ^jk 


♦  6_0  6 


a&  °L (a ,  m)  L(a,  n)[^ 


raLH  ♦ 


*  6a0  *mn  *  jk  , 

as  required.  The  last  step  above  resulted  from  substitution 
of  Equation  (144)  for  aItJ. 

With  regard  to  Equation  (143),  note  that  D"1  reduces 
to  •  form  block  diagonal  in  (“Gj)”1  if  all  the  aIW  happen  to 
vanish.  Thus,  one  may  regard  D”1  as  the  difference  between 
a  "noise”  covariance  matrix  and  a  "residual  bias"  correction 
matrix,  respectively  the  first  and  aecond  terms  of  Equation  (143 


•  (147) 


(148) 
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D.  Some  Expansions  of  the  General  Solution 

The  calculational  efficiency  of  state-vector  estimation 
depends  upon  the  degree  of  difficulty  of  extracting  C-1  in 
the  estimator  Wj  [see  Equation  (70)].  In  "least  squares" 
estimation,  C  is  taken  to  be  diagonal — or  at  least  block 
diagonal,  where  each  block  dimension  is  of  the  order  of  the 
observation-vector  dimension.  Such  a  C  may  be  very  different 
from  D,  if  the  latter  contains  extensive  correlation  terms 
among  residual-bias  errors.  Then  resulting  ephemeris  errors 
will  be  worse  than  the  optimum  (minimum-variance)  errors  when 
C  -  D. 

Of  course,  one  may  choose  C  to  be  a  better  approximation 
to  a  highly  correlated  D,  but  generally  at  significant  cost  in 
calculational  practicality. 

The  purpose  here  is  to  find  explicit,  detailed  represen¬ 
tations  for  Sj,  the  covariance  matrix  of  state-vector 
errors,  for  various  combinations  of  some  possible  structures 
for  C  and  D.  Each  Sj  representation  will  exhibit  a  char¬ 
acteristic  level  of  approximation  to  the  minimum-variance 
form  (when  C  *  D)  and  will  also  afford  a  characteristic 
level  of  calculational  efficiency  not  only  for  Wi,  but 
also  for  Sj. 

Structures  of  D  considered  here  derive  from  the  previous 
section.  Structures  for  C  considered  here  also  derive  from 
the  previous  section,  but  will  occur  in  a  distinctive  nota¬ 
tion  in  order  to  avoid  confusion  with  D. 
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In  the  interest  of  notational  simplicity,  the  ensuing 
expansions  of  Equation  (100)  will  omit  the  subscript  i.  The 
quantities  S,  V ,  and  T  are  still  to  be  understood  as  corre¬ 
sponding  to  the  estimation  epoch  t$.  The  indexing  conven¬ 
tion  of  the  expansions  will  be  the  same  as  in  Section  II1.C, 
i.e.,  it  will  correspond  to  the  observational  epochs  “t^. 

Results  will  occur  as  formulas  for  X*1  and  Z,  where 

(149) 

and 

2  9  (c-1l)tD(c*1T) 

Then,  taking  into  account  the  symmetry  of  C, 

s  -  xzx 


(150) 

(151) 


Also, 

w  -  x(c  1  t)  (152) 

Mote  that  if  C  ■  D,  then  S  •  X.  Thus  X  is  the  covariance 
matrix  of  state-vector  errors  for  the  case  of  optimal, 
minimum-variance  estimation.  Accordingly,  ZX  (or  XZ)  is  the 
matrix  factor  on  X  by  which  S  falls  short  of  optimality. 

Numeric  inversion  of  X-1  to  yield  X  is  not  ordinarily  a 
significant  calculational  problem,  since  X  is  only  a  6  x  6 
matrix  and  need  be  inverted  only  occasionally  [see  discussion 
of  Section  II. E  and  Equation  (102)].  Bence,  one  may  regard  the 
calculational  difficulty  of  finding  w  and  S  as  effectively  as 
that  of  finding  X“ ^  and  Z— the  latter  of  course  con- 

t 

taining  C^T  (see  Equations  (ISO)  and  (152)]. 
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Table  5  gives  seven  combinations  of  structure  assumptions 
for  C  and  D,  along  with  appropriate  analytic  expressions  for 
C-^.  Table  6  gives  corresponding  formulas  for  X”1  and  2. 

The  notes  in  Table  5  pertain,  first,  to  calculational 
ease  of  finding  X”1  and  2,  given  a  set  of  measurement  data; 
and  second,  to  ease  of  recalculating  X-1  and  Z  as  new 
measurement  data  become  available.  Additional  comments  upon 
each  case  will  conclude  this  section. 

Case  1  imposes  no  restrictions  upon  C  or  D,  except  that 
all  error  correlations  vanish  from  one  sensor  to  another. 

These  matrices  are  then  block  diagonal,  respectively,  in 
partitions  °C  and  “d.  The  formulas  for  X-1  and  Z  are  alge¬ 
braically  simple  but  calculationally  complex,  requiring  nu¬ 
meric  inversion  of  the  large  partitions  °Z. 

Case  2  introduces  the  structure  of  Equation  (114)  for  D, 

^  _ 

and  a  similar  structure  for  C  in  terms  of  Ej  (analogous 
to  °fej)  and  (analogous  to  QHj£).  These  structures 

improve  the  physical  interpretabil ity  of  C  and  D,  but  are 

generally  still  not  calculationally  practical  in  that  inversion 

ex 

of  each  of  the  large  partitions  C  is  still  necessary. 

Case  3  reduces  C  to  a  form  block  diagonal  in  QEj,  resulting 
in  a  least-squares  form  of  W.  However,  careful  inspection  of 
the  expression  for  2  in  Case  3  will  reveal  that  the  generality 
of  the  D-structure  places  significant  demands  upon  computer 
memory  in  recalculating  Z  as  new  data  become  available. 


TABLE  6 

STATE  VECTOR  ERROR  MATRICES 
FOR  STRUCTURE  OF  TABLE  5 


Caea  t) 


Case  4  begins  a  different  type  of  structuring,  relative 
to  Case  1  as  a  baseline.  A  pass-multiplet  structure  for  C 
and  D  substitutes  the  problem  of  inverting  the  partitions  aLC 
rather  than  the  cC.  The  aLC  tend,  however,  to  still  be  of 
large  dimension. 

Case  5  is  the  pass-multiplet  analog  of  Case  2  and  is 
calculationally  not  greatly  advantageous  to  Case  2. 

Case  6  introduces  a  structure  for  C  analogous  to  Equa¬ 
tion  (142)  for  D.  The  inverse  of  this  structure  exists  in 
simple  analytic  form  and  makes  calculation  of  C~^  only 
slightly  more  involved  than  for  the  least-squares  Case  3. 
Thus,  the  inversion  Equations  (143)  and  (144)  effectively 
afford  a  limited,  calculationally  efficient  generalization 
of  least  squares  estimation. 

With  Case  6  as  with  Case  3,  however,  the  general  form  of 
D  does  impose  some  calculational  penalties  in  regard  to  Z  as 
new  data  become  available. 

In  Case  7,  both  C  and  D  have  the  structure  of  Equation 
(142).  Calculation  of  Z  is  thereby  considerably  simpler. 
(Further  specializations  of  Case  7  will  appear  in  the 
next  section. ) 

E.  The  Computer  Model  SEEM 

The  purpose  of  this  section  is  to  describe  the  analyti¬ 
cal  basis  of  the  computer  model  SEEM  (see  Chapter  I).  A 
further  purpose  is  to  provide  some  numerical  outputs 
of  SEEM  as  examples  of  calculational  results  obtainable 
using  the  analytical  formalism  of  this  report. 


1 .  Analytical  Basis  of  SEEM 

The  approach  here  is  the  further  specialization  of  Case  7 
of  the  previous  section. 

For  D,  assume  that  each  pass  multiplet  is  a  single  pass. 
Further,  assume  that  a^H  is  the  same  for  every  multiplet. 

Then  one  may  write 


aBnmn 

Djk 


6  e  6 
aB  mn 


jk 


aGm  +  ° 


“) 


(153) 


Make  further  assumptions  about  the  internal  structures 
of  aG™  and  aH.  Assume  that  within  any  given  observation 
the  noise  errors  are  uncorrelated  with  each  other  and  the 
residual  bias  errors  are  uncorrelated  with  each  other.  Then 
both  “Rj1  (see  Equation  (121b)]  and  “Rjk  (see  Equation 
(125)]  are  identity  matrices. 

Then  from  Equations  (119)  and  (122),  the  general  element 
of  QGj  is 


From  Equations  (123)  and  (126), 
aH  is 


the  general  element  of 


(154) 


(155) 
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SEEM  at  present  provides  for  just  one  sensor  type,  radars 
operating  in  altazimuth  coordinates  (see  Section  II. A).*  The 
noise-variance  components  of  Equation  (154)  become 

(v°3)i  '(v°x)2,  11/005  h> 2  ' 

(“<-(X)  2  - 

ori  -o.f  . 

Here  h,  and  p  respectively  denote  azimuth,  elevation,  and 
range,  and  the  symbol  o  on  the  righthand  side  represents  a 
constant  for  a  given  radar.  The  factor  (1/cos  h)2  in 
Equation  (156)  accounts  for  the  loss  in  azimuthal  accuracy 
at  high  elevation  angles. 

SEEM  assumes  the  following  for  the  residual-bias  variances 
of  Equation  (155): 


*  One  may  also,  by  an  input  strategem,  make  SEEM  accommodate 
altazimuth-coordinate  telescopes.  The  strategem  is  to 
assign  a  very  large  value  to  $<^2  of  Equation  (158), 
thereby  assigning  range  measurements  negligible  statistical 
weight. 


(156) 

(157) 
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Ag  ain  the  c  symbols  on  the  righthand  side  are  constants 
for  a  given  radar. 


As  regards  C,  SEEM  provides  two  options, 
a  "least  squares”  option,  is 


The  first, 


mn 


«B/C 
V  “/jk 


The  second,  a  "minimum  variance"  option,  is 

afi/r,  \Inn  _  a£_mn 


‘(CHV) 


jk 


Dj* 


Appropriate  manipulation  of  Case  7  results  of  the  pre¬ 
ceding  section  then  yields  for  the  least  squares  option 


where 


SLS  '  XLs 


(V*  °H  V" 


) 


LS 


and 


3 '  l  1(1  0l?) 


For  the  minimum-variance  option. 


®MV  *  *MV  * 


(162) 


(163) 


(164) 


(165) 


(166) 


(167) 
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where 


and 


Evaluation  of  the  °T?'  in  these  expressions  is  via  Equa- 

J  am 

tions  (98),  since  Tj  there  was  redesignated  at  T^  in  the 
notational  change  of  Subsection  1II.C.1.  On  the  righthand 
side  of  Equation  (98a),  tne  Mj  is  the  "radar"  entry  of 
Table  4.  The  Qj"1  are  those  of  Table  1,  as  specialized  in 
Table  2.  Appendix  B  gives  expressions  for  4>jj  and 
Equation  (98b)  . 

Thus,  SEEM  evaluates  S^g  or  S^v  as  S^  at  an  appropriate 
* 

time  tj.  Equation  (102)  then  yields  S^  at  requisite  pre- 

A  A 

diction  times  t^ ,  k  «  1,  2,  ...  n.  For  each  t*,  SEEM  ex- 
tracts  rS)c  from  S*  (see  defining  Equation  (91  )J  and  then 
obtains  [rSk]lJVW  via  Equation  (103),  evaluating  the  rotation 
matrix  L  as  developed  in  Section  II. A. 

The  diagonal  elements  of  are  the  variances 

°jcU’  °kV*  °kw»  respectively  the  radial,  along-track,  and 
cross-track  component  variances  of  the  ephemeris  vector  at 
tfc.  The  resultant-vector  variance  is 

-  °5k  +  «vk  * 


(166) 


(169) 


(170 


78 


SEEM  outputs  3ojj,  30^,  3jjcv#  and  30^^  for  the  predic- 

A 

tion  times  t^,  k  =  1,  2,  ...»  n. 

As  it  is  now  programmed,  SEEM  does  not  obtain  the  prin¬ 
cipal  axes  and  orientation  of  the  ephemeris  error  ellipsoid 
by  diagonalizing  [rsk]uw$'  etc.  However,  for  near¬ 
circular  orbits  this  ellipsoid  tends  to  be  oriented  along 
the  UVW  coordinates,  primarily  because  of  the  effect  of 
period  uncertainty,  which  makes  the  along-track  error  ordi¬ 
narily  large  compared  to  radial  and  cross-track  errors. 

Thus,  for  near-circular  orbits  one  may  regard  3ojcU, 

3^kV'  and  3°kw  as  approximately  the  principal  dimensions 
of  the  60.8-percent  confidence  error  ellipsoid  (see  Table  A-1). 
The  resultant  30^  is  the  exact  RSS  dimension  of  the  ellip¬ 
soid,  since  the  trace  of  [rsk J  uvw  *s  invar*ant  under  the 
coordinate-frame  rotation  of  diagonalization. 

2 .  Representative  Ephemeris  Error  Results  of  SEEM 

This  subsection  gives  representative  graphical  outputs 
of  SEEM.  All  outputs  given  here  correspond  to  a  satellite 
in  circular  orbit  at  an  altitude  of  400  km.  The  epoch  t  ■  0 
is  that  of  "launch,"  when  the  satellite  initially  appears  in 
orbit. 

Sensors  are  altazimuth  radars,  with  hemispheric  coverage 
down  to  a  minimum  elevation  angle  of  7*.  Maximum  range  is 
set,  not  by  radar  capability,  but  rather  by  horizon  line-of- 
sight  cutoff  at  minimum  elevation.  Sensors  provide  measure¬ 
ments  at  6-second  intervals  while  the  satellite  is  within 
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coverage.  For  "nominal"  measurements,  "noise"  and  "residual 
bias"  standard  deviations  are  those  of  Table  7. 

SEEM  output  graphics  of  Figures  3  through  6  correspond 
to  the  specific  input  conditions  of  Table  8.  Figure  3 
represents  a  pass  through  radar  coverage  about  15  minutes 
after  "launch."  Error  plots  commence  a  few  minutes  after 
the  satellite  exits  coverage,  simulating  a  data-processing 
delay  in  availability  of  estimated  ephemer ides . * 

Note  that  for  this  case,  minimum-variance  estimation 
yields  little  accuracy  improvement  over  least  squares 
estimation. 

The  error  curves  of  Figure  3  exhibit  several  typical 
features  that  deserve  comment.  First,  the  cumulative 
growth  of  along-track  error  is  what  one  would  expect  from 
error  in  estimating  the  satellite  period.  By  contrast, 
the  radial  and  cross-track  errors  repeat  with  each  satel¬ 
lite  revolution.  Thus,  along-track  error  soon  dominates 
the  resultant  error. 

Second,  all  error  components  contain  the  satellite  period 
(100  minutes)  as  a  fundamental  harmonic,  with  minima  at  or 
near  one-period  intervals  from  the  radar  pass.  This  be¬ 
havior  is  not  surprising,  since  the  radar  pass  corresponds 
to  a  point  in  inertial  space  where  position  is  actually 
measured,  i.e.,  where  both  estimated  and  true  orbits  are 
closest  together. 


*  In  Figure  3,  the  width  of  the  "radar  coverage  rectangle" 
denotes  time  in  coverage,  but  the  height  of  the  rectan¬ 
gle  has  no  interpretational  significance. 
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TABLE  7 

“NOMINAL"  STANDARD  DEVIATIONS  OF  MEASUREMENT 


•Coordinate 

Noise 

Residual  Bias 

Azimuth 

®o  -  0.05* 

v  X 

“o  -  0.05* 

SX 

Elevation 

V  «  0.05* 
v  h 

°o  -  0.05* 
e  X 

Range 

°o  ■  50  m. 

v  p 

*0  *  50  in* 

e  P 

TABLE  8 

INPUTS  FOR  SEEM  EXAMPLES 


- — ■  —  ■  - 

Figure 

Std.  Deviations 

of  Measurement 

Pass 

Estimation 

Humber 

Noise 

Residual  Bias 

Spacing 

Method 

3 

(Solid  Curve) 

Nominal 

Nominal 

— 

Least 

Squares 

3 

(Broken  Curve) 

Nominal 

Nominal 

— 

Minimum 

Variance 

4 

Nominal 

Eero 

— 

LS  and  MV 
[Same  Curve] 

5 

Nominal* 

Nominal* 

1/4 

Revolution 

Least 

Squares 

< 

Nominal* 

Nominal* 

1/2 

Revolution 

Least 

Squares 

*  Applies  to  both  radars 


81 


FIGURE  4 

EPHERMIS  ERRORS: 

SINGLE  PASS  WITH  ZERO  RESIDUAL-BIAS  ERRORS 
(Facsimile  of  Computer  Output) 


FIGURE  5 
EPHERMIS  ERROR: 

TWO  PASSES,  QUARTER  REVOLUTION  SPACING 
(Facsimile  of  Computer  Output) 
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Tint  FROM  I.RUMCH  IN  HOURS 


FIGURE  6 
EPHERMIS  ERROR: 

TWO  PASSES.  HALF  REVOLUTION  SPACING 


iwumsM  wim  imi  woiv 

UU1M01MI  Ml  MM3 -0( 
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Third,  the  cross-track  and  radial  errors  clearly  possess 
a  second  harmonic  in  addition  to  the  fundamental.  One  can 
understand  the  structure  of  the  cross-track  error  by  re¬ 
membering  that  the  estimated  and  true  orbital  planes  inter¬ 
sect  on  a  line  running  from  the  (minimum  error)  measurement 
point  through  the  center  of  the  Earth  and  beyond.  One  thus 
expects — and  obtains — a  second  cross-track  error  mini¬ 
mum  where  the  orbits  (nearly)  intersect  180*  away  from  the 
pass  point. 

The  more  complex  structure  of  the  radial  error  probably 
has  to  do  with  the  interplay  of  two  estimated  orbit 
parameters:  the  argument  of  perigee  and  the  orbit  eccen¬ 

tricity. 

Now  compare  Figure  4  to  Figure  3.  "Perfect  calibration" 
of  the  radar  (i.e.,  zero  residual  biases)  yields  a  dramatic 
improvement  in  all  ephemeris  error  components.  One  can 
demonstrate,  in  fact,  that  most  of  the  improvement  results 
from  elimination  of  the  residual  bias  in  elevation. 

Figures  5  and  6  both  contain  a  second  radar  pass,  with 
measurement  errors  the  same  as  for  Figure  3.  The  additional 
data  in  both  cases  lead  to  more  accurate  ephemerides,  as  one 
would  expect. 

Comparing  Figures  5  and  6,  it  is  not  surprising  to  find 
that  the  larger  pass  spacing  yields  a  more  accurate  satellite 
period,  and  hence  reduced  along-track  figure  error.  The  ra¬ 
dial  error  is  also  slightly  better  for  diametrically  opposed 
points  on  the  orbit. 
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On  the  other  hand,  the  quarter-orbit  spacing  of  the  passes 
should  stabilize  the  estimated  orbital  plane  (which  must  con¬ 
tain  the  center  of  the  Earth)  better  than  half-orbit  spacing. 
Thus,  one  may  understand  the  reduced  cross-track  error  of 
Figure  5  relative  to  that  of  Figure  6. 

In  all  of  the  foregoing  examples,  the  pass  geometry  is 
such  that  the  satellite  flies  almost  directly  over  the  radar. 
Reference  5  contains  SEEM  output  examples  for  other  pass 
geometries,  for  lower  minimum-elevation  angles,  and  for  as 
many  as  three  passes. 
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APPENDIX  A 


COVARIANCE  MATRICES:  DEFINITION  AND  PROPERTIES 

This  Appendix  provides  a  review  of  the  definition  and 
primary  properties  of  covariance  matrices.  The  first  sub- 
-section  gives  properties  not  restricted  to  any  particular 
probability  density  function.  The  second  subsection  treats 
covariance  matrices  for  the  special  case  of  normal  (Gaussian) 
distributions. 

1 •  Definition  and  General  Properties 

Consider  an  arbitrary  random  n-vector 


ft*  »7 


(A-l  ) 


where  fi  is  the  true  value  of  w  and  rj  is  some  random  error, 
such  as  measurement  error. 

Let  the  probability  density  function  of  ri  be  p{>7).  The 
probability  density  of  the  ith  component  is  then 


pi(V  “  /*••/  PW  *  (dnk>  . 


(A-2) 


(The  subscript  on  p^  indicates  that  its  functional  form  may 
depend  upon  the  value  of  i.)  The  joint  probability  density 
of  ni  and  nj  (where  i/j)  is 


pij(T1i'V 


*  /**•/  p(i?)  n  (dnt)  . 

—  —  k?i  * 


(A-3) 


PNEGBKW  MS  BJUK-MOT  Flli® 


Clearly  this  function  is  symmetric: 


Pij(VV  *  Pji<VV  * 

With  these  definitions,  the  probability  that  n*  lies  in 
the  interval  (a,b]  is  then 

b 

l  Pi<Vd”i  • 

The  probability  that  hi  lies  in  la.b],  while  also  n  lies 
in  tc,d]  is  3 

d  b 

/  J  . 

Now  define  the  expectation  value  of  any  function 
f  (17)  as 


Then 


£{f  (17) ) 


n 

«(»?)p(i7)  n  (dnv)  . 

x- 1  * 


“_/  nip(ni)dni 


*  , 


the  mean  of  n*.  ‘  If  fjm  0,  then  17  is  said  to  be  unbiased. 
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(A-4) 


<A-5) 


<A-6) 


Also, 


e{  (ni-“i)2) 


(T1i*^i)2pj  (T)i)drii 


i 


the  vat^anct  of  m«  Further, 


(A- 7 } 


'-i  (nj’~j)Pij  (ni,n^)dr,idrij 


s  p .  .0.0. 

i]  i  3 


(A-8) 


the  cova**ance  of  ni  and  nj*  Define  quantity  pij  as  the 
coKKilation  coiliiint  of  n i  and  ni.  A  theorem  exists  that 


-1  ^  pij  i  • 


(A-9 ) 


Finally,  define  the  covatiance  matKix  Cn  a  matrix  with 
diagonal  elements 


<cVu  *  *i 


(A-lOa) 


and  off-diagonal  elements 


<cn’ij  ’  *ij Vj  . 


(A-lOb) 


That  ia. 


C  -  m-H)')  . 


(A-U) 
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Note  that  one  may  estimate  C from  a  set  of  data  sam¬ 
ples  (I7)fc»  k  ■  1,  2,  . ..,  K  by  replacing  the  integrals  of 
Equations  (A-6)  and  (A-7)  by  sums  over  the  index  k,  providing 
that  the  joint  probability  functions  P^fni.nj)  sre  known 
(or  assumed).  (This  report  does  not  pursue  further  the  sub¬ 
ject  of  sampling  and  subsequent  estimation  of  C  . ) 

h 

Clearly  is  symmetric.  If  all  the  r>i  are  linearly 
independent,  then 


and  one  can  show  that  C  is  positive  definite — i.e.,  its 

-i  '' 

inverse  Cn  •  exists. 

Now  consider  the  transformation 

tm  F17  , 

where  (  is  an  m-vector  (m<n)  and  r  is  any  m  »  n  matrix  that  is 
not  a  function  of  17. 

A  functional  dependence  of  f  on  fj ,  such  as  that  in 
Equation  (A-13),  means  that  for  every  point  occurring  in  ij- 
•pace  a  corresponding  point  exists  in  ( -space.  Thus,  for 
any  arbitrary  volume  in  17-space,  the  same  number  of 
points  must  occur  in  the  corresponding  volume  in 
{-space: 


***drin 


-  /•••/pn(i7)dh1 


( 


Here  the  subscripts  on  the  probability  densities  denote  that 

their  functional  forms  are  in  general  different.  But  since 

V  is  arbitrary  in  Equation  (A-14),  it  follows  that 

n 

Pc  (f)d;1»««d;|n  «  p^  (Wdri^'driH  . 

This  is  a  fundamental  invariance  relation  for  functional 
transformations  among  random-vector  spaces. 

It  now  follows  that 


(A- 1 5 ) 


that  is. 


Similarly, 


E{£ )  5  /•••/  f  P, (£)d;  **.dt 

*00  mod  V  X  X 


"P"  T'*' 

=  /•••/  Trj-p  (rj)dn,»‘-dT)n 
_cc  _ce  n  A  n 


f*e  {17}  , 


1  m  FT}  . 


CC  “  F*VF  • 


(A- 1 6 ) 


(A-1 7 ) 


(A-18) 


Note  that  neither  Equation  (A-17)  nor  (A-18)  is  necessarily 
valid  if  P  ■  F ( 17 ) . 

According  to  Equation  (A-17),  if  17  is  unbiased  then  ( 
will  also  be  unbiased.  Equation  (A-18)  tells  how  to  calcu¬ 
late  the  covariance  aiatrix  in  (-space.  One  can  show  that 


1 


the  symmetry  property  of  is  passed  on  to  C^.  Further¬ 
more#  if  exists,  and  if  F  is  of  rank  m,  then 

exists  and  is  also  symmetric. 

Consider  the  special  case  in  which  the  transformation 
matrix  is  Q,  where  Q  is  an  n  x  m  matrix  of  full  rank  such  that 


that  is. 


QQ  -  I  , 


Q+  *  C"1  . 


(A-1 9a ) 


(A-1 9b) 


Examples  of  such  matrices  are  those  that  perform  orthogonal 
rotations  upon  the  space  17.  Of  particular  importance  is  the 
rotation  Q  which  yields  a  diagonal  fdiagonalization  of 

cn->. 

Finally,  one  can  show  that  when  F  *  Q,  the  determinants 
of  Cg  and  Cf,  are  equal: 


|CCI  -  ICJ  . 


(A-20 ) 


Moreover,  the  trace  of  C ^  is  invariant  under  the  Q- 
trans format ion: 


(A-21) 


Thus,  if  a  variance  vector  e  is  defined  with  components  oy, 
®2»  ••••  ®n»  then  Equation  (A-20)  states  that  the  length 
of  0  is  invariant  under  coordinate-system  rotations. 


2.  Properties  When  the  Density  Function  Is  Normal 

Properties  of  17  end  treated  thus  far  do  not  depend 
upon  any  particular  assumptions  about  the  form  of  p ( 17 ) . 
Consider  now,  however,  the  following  prbblem. 

Referring  to  Equation  (A-1),  suppose  a  measurement  of fi 
is  attempted,  yielding  the  measurement  »  because  of  noise 
corruption  17.  Suppose  one  knows  both  ?)  and  C^.  Let  the 
estimated  value  of  fi  be 


Then  the  mean  of  a  large  number  of  measurements  will  be 


(A-22) 


P*  «  V 


(A-23) 


The  problem  is  now: 

a.  What  is  the  probability  P  that  f/*  will  fall 
within  some  prescribed  volume  V  about  the 
point  /#? 

b.  What  is  an  appropriate  prescription  for  V? 
As  regards  Question  a.,  obviously 


p  «  /•••/p(i7)dn  •••dn_  . 

V  in  (A-24) 

Evaluation  of  P  requires  knowledge  of  the  integrand  within 
V.  However,  17  and  C^,  the  first  and  second  moments  of 
p(»7),  do  not  in  general  completely  specify  p(97) — which  may 
possess  higher,  independent  moments.  Hence,  additional 
assumptions  about  p(Tj)  are  necessary  in  order  to  evaluate  P. 


A- 9 


Moreover,  one  can  expect  that  the  answer  to  Question  b. 
will  be  sensitive  to  the  assumed  form  of  p (»)). 

Assume  now  that  p(i))  has  normal  (i.e.#  multivariate 
Gaussian)  form: 


j  tcn'1  (n-n> 

pW’ '  ‘ 

This  probability-density  form  has  widespread  utility.  Since 
the  form  is  completely  specified  by  and  ,  one  can  expect 
a  mathematically  fruitful  investigation  of  the  stated 
problem. 

The  investigation  proceeds  by  considering  the  behavior 
of  ( 77 )  under  coordinate-system  rotations,  reverting  now 
to  the  notation  of  Equation  (A-15).  Now  under  a  coordinate 
rotation 

*m°1  (A-26) 

the  Jacobian  of  the  transformation  is  the  absolute  value  of 
the  determinant  |Q|,  and  ao 

d;l**’dtn  *  li0lidrV”dnn  *  (A-27 ) 

Then  by  Equation  (A-15) 


P«)  *  “7.  PW 


II  011 

-  J(»»-?)+Cn"1(»7-W) 


(2i0n/2|jQlNO1/2* 


(A-28 ) 


(A-29 ) 


! 
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Then,  using  Equation  (A-20), 


Il0||-|cnl1/2  -  |c;|l/2. 

Also,  using  the  rotation-matrix  property  of  Equation  (19a), 

-  (H-7) +oo+cli*1oQt(i) -W) 

-  «-«)  +cc_1  «-T>  • 

Hence 


P;(f) 


a  *  v 


(A-32) 


which  is  of  the  same  functional  form  as  Equation  (A-29). 

One  can  also  readily  infer  that  under  coordinate-system 
translations,  the  form  of  Equation  (A-29)  is  also  invariant. 
Let  this  form  be  denoted  as  p(^,Cn;n).  Thus,  the  mean  and 
the  covariance  matrix  play  the  same  roles  in  the  probability 
density  function,  regardless  of  coordinate-system  selection. 

Now  assume  that  Q  is  chosen  such  that  is  diagonal. 

Then 


p({,C;;f) 


(2x)n/2(o  *o 

12  n  ( 


•  a 


u.-?, )2  (;2-T2)2  un-;n>2  ' 

PA  <°a\  <°n\  J  (M31) 


A-ll 


n 

n 

i=l 


(A-33b) 


(oi)c 


i  “i-h>2 

2 — 

s  <°A 


(i.e. ,  the 
covariance 
density  is 
components 


factor  1/  ^occurs  n  times).  Thus,  if  the 
matrix  is  diagonal,  the  multivariate  probability 
the  product  of  univariate  densities — i.e.,  the 
of  f  are  statistically  independent. 


This  result  does  not  hold  in  general  for  non-normal 
forms  of  p.  However,  one  can  show  that  statistical  inde¬ 
pendence  always  implies  a  diagonal  covariance  matrix. 


This  completes  the  groundwork  appropriate  to  address 
Question  b. 


Suppose,  in  Equation  (A-25),  the  variable  expression  in 
the  exponent  is  set  equal  to  a  constant  a2,  where  a  >  0: 


»  a2  .  (A-34) 

This  is  the  equation  of  an  ellipsoid  in  the  n-fold  hyperspace 
of  ri.  For  example,  if  n  «  2,  Equation  (A-34)  becomes 


(r^-n)2 


-  2p 


(tJj-tTj)  <h2-rT2)2 


12 


a  0 
12 


.2U-P12J) 


(A-35 ) 


Suppose  the  integration  volume  V  is  taken  to  be  this  n- 
dimensionai  "error  ellipsoid.”  Then  the  questions  are:  What 
are  the  properties  of  this  error  ellipsoid,  and  to  what  ex¬ 
tent  is  it  described  by 


A-12 


The  center  of  the  ellipsoid  in  tj- space  is  obviously  at 
If.  By  Equations  (A-1 )  and  (A-22), 


(A-36) 


so  that  in  space  the  ellipsoid  center  is  at  fl ,  which  is 
certainly  desirable. 

Prom  Equation  (A-34)  one  can  see  that  the  surface  of 
the  error  ellipsoid  is  a  surface  of  constant  probability. 
For  a  rotation  of  the  coordinate  system  to  a  new  t-set 
[see  Equation  (A-31)],  the  form  of  Equation  (A-34)  is 
clearly  invariant.  The  relationship  of  the  ellipsoid  to 
the  new  probability  density  [Equation  (A-32))  is  also  in¬ 
variant.  Thus,  taking  Equation  (A-20)  into  account,  one 
can  see  that  the  numerical  value  of  the  probability  density 
on  the  ellipsoid  surface  is  invariant  under  rotations. 

Since  a  is  arbitrary,  it  follows  that  P — given  by  Equation 
(A-24) — is  also  invariant  with  respect  to  coordinate-system 
choice.  Thus,  the  error  ellipsoid  can  meaningfully  serve 
as  a  "confidence  volume." 


Assume  now  that  the  rotation  Q  has  been  selected  so  as 
to  make  diagonal.  Then,  from  Equation  (A-33a) 


—  2  —  2 
(t  -t  >  U  -C  ) 

11^  2  2 

4  - a -  ♦  ••  •  + 


“A 


‘v’V 


«SrV 
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The  principal  half-axes  of  the  ellipsoid  are  clearly  a(ci)^» 

«(o2)£»  ...»  alo n>£»  The  ellipsoid  orientation  in 
Tj-coordinates  is  specified  by  the  requisite  diagonali- 
zation  matrix  Q.  That  is,  Q  is  the  set  of  direction  cosines 
between  the  17-axes  and  the  (-axes  along  which  the  ellipsoid 
is  oriented. 

Thus,  C  (and  of  course  also  C-1 )  contains  information 
that  specifies  both  the  shape  and  orientation  of  the  error 
ellipsoid  (0  is  derivable  from  Cn ) .  The  selectable  parameter 
a  is  the  scaling  factor  of  the  principal  dimensions  of  the 
ellipsoid.  Thus 

V  «  V(Cn,a,n)  .  (A-38) 

The  "1c  -  confidence  volume"  is  then  V(C^,1,n),  the  "2o  - 

confidence  volume"  is  V(C  ,2,n),  etc. 

h 

One  can  extract  a  certain  amount  of  information  from  Cn 
and  Cn~1  without  diagonal ization.  The  standard  deviation 
of  the  resultant  of  17  is  the  trace  of  C^: 


(see  Equation  (A-21)).  The  are  the  reciprocals  of 

the  intercepts  of  the  1o  -  ellipsoid  with  the  coordinate  axes 
of  the  U-system. 

Question  a.  can  now  be  answered  by  integrating  Equation 
(A-24)  in  the  most  convenient  coordinate  frame,  i.e.,  the 
diagonal  frame  (  for  (  ■  0: 
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*■/•••/  P(o,c.,£)d;  •••« 

V(C;,a,n)  4  1  » 

n  — —  x  •  ^ 

*/••*/  H  — — e  2  *  dx.  , 

Sphere  of  i«l  /5rT  1 

Radius  a 

using  Equation  (A-33b)  with  a  change  of  coordinates 


(A-40) 


X  =  — i —  .  (A-41) 

‘#l>t 

Thus  P  for  n  *  1 ,  n  *  2,  n  ■  3,  etc.,  is  respectively  the  proba¬ 
bility  associated  with  a  linear  confidence  interval  of  length  2ao, 
a  circular  confidence  area  of  radius  ac,  a  spherical  confidence 
volume  of  radius  ao,  etc.  Table  A-l  gives  numerical  values  of 
P  for  some  commonly  occurring  values  of  n  and  a.  This  completes 
the  investigation  of  Question  a.  for  n-vectors. 

Suppose,  now,  that  Questions  a.  and  b.  are  extended  in 
scope  to  read  as  follows: 

(1)  What  is  the  probability  Pro  that  m  specified  com¬ 
ponents  of  ft*  will  fall  within  some  prescribed 
volume  Vjj,  about  (i.e.,  referring  to  the 
specified  m-dimens ional  subspace  of  ft  ,  m  <  n)? 

(2)  what  is  an  appropriate  prescription  for  Vm? 

As  before,  these  questions  will  be  considered  only  for  the 
case  of  normal  probability  density  p (17,0^? q). 
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TABLE  A-1 

SOME  USEFUL  P  VALUES 


To  simplify  the  ensuing  notation,  assume  for  the  moment 
that  the  first  m  components  of  f/*  are  those  specified.  Then 
the  probability  density  for  the  m  components  is 


P(0  *  /**•/  P<^Cn»^drWl‘*‘drin  ' 


(A-42) 


pm  *  /;/•/  p(VdV**dni 


(A-43) 


The  solution  of  the  subspace  problem  then  hinges  upon  a 
fundamental  theorem  of  multivariate  normal  probability.  One 
can  prove  that 


P  (*?_)  *  - =773 — 

™  (2r)In/2|C 


i _ _  -»L)  +c  ~1  w-t)  ) 

©  »  in  m  m  m 


•  p(VtV) 

Here  the  covariance  matrix  C  .  is  the  subset  array  of 

n  in 

corresponding  to  the  m  specified  components  of  ft*. 

Since  the  forms  of  Equations  (A-42)  and  (A-43)  are  iden¬ 
tical  to  those  of  Equations  (A-25)  and  (A-24),  respectively, 
all  of  the  previous  results  for  V  and  P  of  the  full  n- 
vector  immediately  follow  where  one  substitutes  rim  and  C 

nm 

for  V  and  C 


(A-44) 


Extension  to  any  m  components  of  ft*  is  obvious,  since 
the  numbering  of  the  components  is  arbitrary.  One  must,  of 
course,  select  corresponding  components  in  forming  f}m  and 

cnm* 

Finally,  up  to  this  point  the  analysis  has  assumed  that 
both  Cn  and  7}  are  known.  If  only  rj  is  known,  however,  one 
can  find  ft*  but  not  P.  If  only  is  known,  one  can  find  P 

but  not  n*  •  These  comments  also  apply  to  m-fold  subspaces 
of  these  n-vectors. 
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PROPAGATION  ERROR  MATRIX 

The  first  of  the  following  subsections  states  the  problem 
of  determining  a  matrix  representing  the  error  in  the  propa¬ 
gation  matrix  [Equation  (B-2)].  The  second  subsection  gives  the 
solution  approach  and  the  solution  itself.  The  final  sub¬ 
section  contains  intermediate  steps  of  the  analysis. 

1 .  The  Propagation  Error  Matrix  Problem 

The  objective  here  is  to  find  the  effect  of  error  in  the 
state  vector 


(B-l  ) 


upon  the  propagation  matrix 


f  0  0  ]  g  0  0 

0  f  0  ,  0  g  0 

0  0  f  |  0  0  g 

f  o  0  ;  g  0  0 

0  f  0  ,  0  g  0 

0  0  f  J  O'  0  g 


(B-2) 


Here  (see  Reference  B-1)* 


*  -  1-  (X-cos  i) , 


(B-3a ) 


*  Equations  (B-3a  to  3d)  are  among  Equations  (7P18)  and  (7P19) 
with  the  convenient  abbreviations  p  =  cD,  q  s  d0//T. 
Equations  (B-4)  to  (B-6)  are  among  Equations  (7B5)  -  (7B9). 
Equation  (B-8)  is  Equation  (7P14),  and  Equation  (B-9)  is 
obtained  by  integrating  Equation  (4G4)  with  E  s  E-Eo  end 
n  defined  by  Equation  (4G7),  identical  with  Equation  (B-7) 
here. 
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with 


and 


g  «Q|jsin  I  ♦  •&-  (1-cos  E) ; 

(B-3b) 

f  -  -”(^)tln  =' 

(B-3c ) 

g  *  1-  (l-cos  1) > 

(B-3d) 

-(n  -¥f. 

(B-4) 

p  -  l-  i  i 

(B-5 ) 

>o'fo  , 

q'-75T’ 

(B-6 ) 

n  -  /y  a'5''2. 

(B-7 ) 

r  *  a(l  -  p  cos  E  ♦  qsinE), 

(B-8) 

1  rft 


•  -jj- (E  -  P*inE  ♦  q (1-cos  EJJ. 
For  Earth  satellites /jT s  3.78808  *  10*  km3/,2/min. 


(B-9) 


B-4 


Specifically,  the  problem  is  to  obtain  an  explicit  repre 
sentation  for  *tt0»  the  Propagation  truor  matnixi  as  defined 
by 


64>^ 

tt 


Yaa 

tt 


6  * 


2.  Solution  Summary  and  Results 

Using  Equations  ( B— 1 )  and  (B-2),  the  lefthand  side  of 
Equation  (B-10)  is 


r  6f  + 
o 

r  + 
o 


where 


3  f  _  f  3  f  3f  3  f 

TT~  ~  3x_  7y~  Tz~  ' 
o  o  o  o 

~  # 


the  quantities  xQ,  y0,  zQ  being  the  components  of  r0;  etc. 
Comparing  the  Equations  (B-11)  and  (B-10), 


( B — 1 0 ) 


»(B- 1  1  ) 


(B- 1 2 ) 


(B-13) 
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The  next  subsection  shows  that 


where 


df 

fllro 

+  f.2r  +  , 

df 

t  ,  •  t 

d  f 

f21o 

22  o  1 

o 

II 

H 

(?) 

in  +(t!) 

o 

1  1  1 

i  df 

f12  - 

It 

'  na  / 

1  7q  ' 

f  21  * 

(=*! 

1 3f 

1  Tq  ' 

/  2 

\  »t  J  Jr°\ 

f  22 

to 

) 77  \7?) 

df 


■(5) 


df 


(B- 1 4a  ) 


(B-Hb) 


(B- 1 5a  ) 


(B-15b) 


(B- 1 5c  ) 


(B-1 5d  ) 


Identical  forms  exist  in  g,  f,  and  g.  Using  Equations  (B-14) 
and  their  analogs,  the  required  result  is 


‘fllVro 

4  f12 

,  •  t 

+JnVro_ 

4  *12 

f . .  ff  •  r  * 
11  o  o 

•  •  *  f 

• 

u  ®11  ro*  ro 

4  *12 

:  t 


f2lVro  4  f22V 


:  4«21r‘o‘rot  4  *22  V 


£21ro*ro+  4  ^22  ro* 


Vro  .4«2lVro 


«22ro* 


(B-16) 
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The  partial  derivatives  of  Equations  (B-15)  and  of  its 
analogs  in  g,  f,  and  g  are  given  in  Table  B-1,  wherein: 


K  5  (f  •  3a)n(l  '  V 

Ep  =  -p  sin  E  ,  (B-1 7b ) 

Eq  =  “  -f-(l“COS  E);  (B-1 7c) 

w  =  p  cin  £  +  q  cos  E  .  (B-1 8) 


Note  that  expressions  for  the  derivatives  of  f,  g,  f, 
and  g,  with  respect  to  each  of  the  components  of  r0  and 
rc,  are  available  in  Reference  B-2.  There,  development  of 
derivatives  proceeds  from  expressions  of  the  general  form 
of  Equations  (B-14)  and  (B-15)  here,  except  that  the 
intermediate  variables  D0,r0,  and  1/a  are  used 
rather  than  p,  q,  and  a.  One  can  prove  that  the  results 
here  are  identical  to  those  of  Reference  B-2,  after  cor¬ 
recting  certain  misprints  in  the  latter.* 

3.  Derivation  of  Partial  Derivatives 

This  subsection  contains  derivations  of  Equations  (B-14), 
(B-15),  and  results  of  Table  B-1  with  defining  Equations 
(B-17)  and  (B-18). 

_  8M 

*  In  Reference  B-2,  Equation  (15C20),  the  term  THT  has  been 
omitted  [see  Equation  (15C5)].  In  Equation  (15C53)  for 
gr,  the  first  denominator  should  read  rDr^  rather  than 

r©3. 
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TABLE  B 1 

PARTIAL  DERIVATIVES 


Referring  to  the  first  footnote  of  this  appendix. 
Reference  B-1  shows  that 


p  «  e  cos  E_  , 
o 


(B-1 9a ) 


q  ■  e  sin 

o 


(B-1 9b) 


Here  e  is  the  orbit  eccentricity  and  E0  is  the  eccentric 
anomaly  at  epoch  t0.  One  will  find  that  of  the  para¬ 
meters  a,  p,  q,  n,  and  r,  only  three  can  be  independently 
chosen. 

Let  these  three  be  a,  p,  and  q.  Then  if  £  is  any  position 
or  velocity  component. 


(B-20 ) 


with  analogous  relations  in  g,  f,  and  g.  From  Equations 
(B-4)  through  (B-6)  one  can  determine  that 


^  x°  ’ 


(B-21a) 


o  r 


■(l+p)x0  , 


(B-21b ) 


xo  * 


(B-21c) 
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with  analogs  in  y0  and  z0.  Similarly 


3a 


i£. 


« 


(B-22a ) 


(B-22b ) 


(B-22c) 


with  analogs  in  y0  and  z0.  Equations  ( B— 14)  and  (B-15), 
and  analogs  in  g,  f,  and  g,  follow  immediately  by  combining 
Equations  (B-20),  (B-21),  and  (B -22). 

Next,  note  that 


3n  _  3n 

Ja  ”  2a  '  (B-23a ) 


3n 


0 


(B-23b ) 


3n 


0  . 


Also,  for  an  arbitrary  variable  n, 


(B-23c) 


(B-24) 


B-10 


by  Equation  (B-5).  Hence 


8 

77 


ft)  ■  *  ■ 
ft)  ■  fey 


3 


3  1  \  -  0  . 


)q\r 


From  Equations  (B-8)  and  (B-18), 


&(■£■)  '  *(-r)  w  H 

^(■f)  ■  '(-r)  (-cot  1  4  w  H)  > 


3 

Sq  r 


T  "  '(-f)  (,in  1  4  w  H)  • 


Differentiation  of  Equations  (B-9)  yields  Equations  (B-17), 
where 


(B-25a ) 

(B-25b ) 

(B-25c) 

(B-26a  ) 

(B-26b ) 

(B-26c ) 


SE  . 

71  “ 


(B-27a) 


B-ll 


can 

and 


it  - 
*P  S 


( 


it  - 

*q  - 


( 


With  the  aid  of  Equations  (B-23),  (B-25),  and  (B-26)  one 
then  differentiate  Equations  (B-3)  with  respect  to  a,  p, 
q  to  obtain  the  results  in  Table  B-1. 


B-27b) 


B-27c ) 


B-l  2 


REFERENCES 


B-1.  Herrick,  S.,  A-4  tKodynamicA ,  Vo  l. 
London,  1971. 

B-2.  Herrick,  S.,  Att* odynanici ,  Vo t. 
London,  1972,  pp.  51-55. 


1,  Van  Nostrand  Re  inhold, 

2,  Van  Nostrand  Re  inhold. 


B-13 


APPENDIX  C 

INVERSION  OF  A  MATRIX 


THE  PREDICTION 
OF 

SATELLITE  EPHEMERIS  ERRORS 
AS  THEY  RESULT  FROM 
SURVEILLANCE- SYSTEM  MEASUREMENT  ERRORS 


mamm  tv*  ■a*-**0*  nu* 


APPENDIX  C 


by  standard  manipulative  techniques.  Cramer's  Rule  then 
holds  for  finding  the  inverse  matrix,  whose  corresponding 
partitions  will  again  be  members  of  the  field. 

Now,  the  set  of  n  *  n  invertible,  diagonal  matrices  plug 
the  null  matrix,  forms  a  field  under  matrix  addition  and  mult 
plication.  Both  Gj  and  H  are  members  of  this  field,  assuming 
for  the  moment  that  H  as  well  as  Gj  is  invertible.  Thus 

|C|  •  -  |MIa  l  I,J  -  1# 2, •  •  •  ,n 

an  n  x  n  matrix  equation  wherein  all  quantities  are  to  be  ex¬ 
pressed  in  terms  of  Gj  and  H.  Here  |c |  is  the  determinant 
of  C.  The  quantity  |aJij  is  the  cofactor  of  Cjj,  and  (C’Mjj 
is  the  general  element  of  C"1. 

One  may  us<>  a  recursive  approach  to  evaluate  the  deter¬ 
minant  |C|.  Let 

CN  *  C  • 

Then 

Cx  -  CA  +  H  , 

and  for  N>1  let  CN-<|  be  the  matrix  of  the  first  N-l  rows  and 
columns  of  C|j. 

Now  one  may  expand  |cN j  via  the  last  row,  and  after  some 
rearrangement  of  rows  and  columns  obtain  the  following: 


S  *  h,-Icn-ii 


H  H  H 

% 

•  •  • 

•  •  • 

•  •  • 

G^+H  H  H 

H  H 

H  H  H 

H  H  H 

•  •  • 

•  •  • 


S-34*  H  H 
H  H 


H  H  B 

H  H  H 

•  •  • 

•  •  • 

•  •  • 

G^3«H  H  H 


H  G^_2*W  H 


One  may  evaluate  the  determinants  here  by  subtracting  the 
last  row  from  each  of  the  other  rows  and  then  expanding  via 
the  last  column.  Equation  (C-6)  then  becomes 


5/N-l  \N-1  , 

IS’1  *  <S  +  GI  • 


By  evaluating  { C 2  I  and  then  | C 3  | ,  one  may  infer  that  |C^ 
probably  has  the  form 

lc»'  -  X1  *  \Wl)  • 


(C-6  ) 


One  may  easily  prove  this  result,  for  N>1 ,  by  substitution 
into  Equation  (C-7). 

Clearly  fc^l  is  n  *  n  diagonal,  and 

lc»'u  ■  +  "«JL(i/(c*  U  ]  • 


By  Equations  (C-2)  and  (C-3),  none  of  the  |CN|jj  vanishes, 
»o  that  |CnM  exists.  Hence  by  Equation  (C-4)  all  of  the 
Cjj”1  exist  and  so  C" 1  exists. 

Now  in  Qeneral  |du|  is  an  (N-l)  x  (N-1)  matrix. 

Define 

K.l  -  Ur.il  . 


(C-9  ) 


(C-10) 


where  for  consistency  of  the  preceding  results 

*  1  •  (C-11 ) 

o 

the  n  *  n  identity  matrix.  The  values  of  MijL#  corresponding 
to  N-2,  are  clearly 


|Aix|  -  G2  4  H  I  |A12I  *  -H  ; 

1^21 1  x  “  "H  *  U22l  *  Gi  +  H 


For  N>2,  one  can  find  l^ula-i  by  again  using  a  recursion 
approach.  For  the  cases  1* J,  one  may  find  via 

exactly  the  approach  used  previously  for  |Cfl|: 


'*«'«  *  *  "X6 k*1  -  ^r1)  • 


(C-T2) 


Observe  that  this  result  holds  for  N>1. 

For  1/J ,  a  typical  | j  -j  for  large  N  is 


Gj4H 


^N-2,N-3*n-1  “  l-1> 


N-2(„i)N-3 


H 

H 

H 


...  gn.34H  h  h 
...  H  H  H 


H  H  Gn4H 


(C-13) 


C-7 


By  interchanging  the  last  two  rows  and  then  the  last  two 
columns,  one  can  bring  the  determinant  here  into  the  form 
of  the  determinants  of  Equation  (C-6)  and  use  the  same 
evaluation  procedure  as  before.  One  can  then  see  that 


(' 


*1J 


-B 


(C-U  ) 


valid  for  N>2. 

One  can  then  easily  verify  that 


(015) 


where  cjj  is  the  Kronecker  delta.  This  holds  for  N> 1 , 

1*1,..., N  and  J*1 , . . . ,N.  This  equation  also  holds  even  when 
B  is  not  invertible,  since  that  particular  property  of  B — assumed 
earlier — was  not  actually  used  in  the  derivation.  Finally, 
one  can  decompose  this  diagonal-matrix  equation  into  a  set  of 
scalar  equations  by  inspection,  with  the  aid  of  the  relation 
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